Modelling the Effect of Communicative Attempt (H1) and Answer Similarity (H2) on Effort
2024-10-17
Introduction
In this script, we will be modelling the causal relation between effort and correction to confirm/reject our hypothesis.
These are:
H1: In corrections, people enhance some effort-related kinematic and/or acoustic features of their behaviour relative to the baseline.
H2: The enhancement depends on similarity of the guesser’s answer and the original meaning. More similar answer will require/result in smaller enhancement (but still enhancement) than less similar answer.
To assess the design and performance of the model, we will use synthetic data that we create to have certain interdependencies and where the effects have pre-defined values. We use this approach instead of using our dyad 0 data, because the pilot data do not include enough data to have a sensible testing sample.
Since we cannot predict the outcomes of the XGBoost modelling (in the previous script) and therefore do not know the variables we will model after preregistration, our model only assumes a continuous variable but it is otherwise free of any other assumptions.
Setting up the environment
library(here)
library(dplyr) # for data-wrangling
library(lme4) # for linear mixed-effects models
library(tidyr) # for reshaping data (if needed)
library(ggplot2)
library(tibble)
library(rcartocolor)
library(ggdag) # for dag
library(dagitty)
library(ggplot2) # bayesian stuff
library(patchwork)
library(bayesplot)
library(brms)
library(beepr)
library(bayestestR)
library(tidyverse)
# current folder (first go to session -> set working directory -> to source file location)
parentfolder <- dirname(getwd())
datasets <- paste0(parentfolder, '/09_Analysis_Modeling/datasets/')
models <- paste0(parentfolder, '/09_Analysis_Modeling/models/')
plots <- paste0(parentfolder, '/09_Analysis_Modeling/plots/')Loading in our data
#TODOBecause current pilot data are quite likely unsifficient in terms of power, and would hinder testing the design of the statistical models, we will create synthetic data with assumed coefficients for all predictors we plan to add to the model.
Creating synthetic data
# Set seed for reproducibility
set.seed(0209)
# Set coefficients
b_exp_vocal <- 0.6 # Vocal has lower expressibility
b_exp_multi <- 1.5 # Multimodal has higher expressibility
b_bif <- 1.15 # More extroverted → more effort
b_fam <- 1.10 # More familiarity → more effort
b_exp <- 1.20 # More expressibility → more effort
b_multi <- 0.70 # Multimodal → slightly reduced effort
b_comatt2 <- 1.50 # Effort increases for second attempt
b_comatt3 <- 0.50 # Effort decreases for third attempt
b_prevan <- 0.50 # Higher previous answer similarity → less effort
# Define participants and concepts
n_participants <- 120
n_total_concepts <- 21 # Each participant works with all 21 concepts
n_modalities <- 3 # Gesture, vocal, combined
# Generate participant-level data
participants <- 1:n_participants
Big5 <- runif(n_participants, min = 0, max = 2)
Familiarity <- runif(n_participants, min = 0, max = 2)
# Generate expressibility scores for each concept across modalities
expressibility_matrix <- matrix(runif(n_total_concepts * n_modalities, min = 0, max = 1),
nrow = n_total_concepts, ncol = n_modalities)
# Initialize data storage
final_data_synt <- data.frame()
# Define function to simulate participant data
simulate_participant <- function(participant_id) {
participant_data <- data.frame()
trial_number <- 1
for (concept_id in 1:n_total_concepts) {
# Assign a random modality
modality <- sample(c("gesture", "vocal", "combined"), 1)
# Get expressibility score based on modality
expressibility_score <- ifelse(modality == "vocal", expressibility_matrix[concept_id, 1] * b_exp_vocal,
ifelse(modality == "gesture", expressibility_matrix[concept_id, 2],
expressibility_matrix[concept_id, 3] * b_exp_multi))
# Base effort before adjustments
base_effort <- b_bif * Big5[participant_id] +
b_fam * Familiarity[participant_id] +
b_exp * expressibility_score +
rnorm(1, mean = 1, sd = 0.5)
# Adjust effort for multimodal condition
if (modality == "combined") {
base_effort <- base_effort * b_multi
}
# Sample the number of communicative attempts (CommAtt)
adjusted_prob <- c(1 - Familiarity[participant_id],
1 - Familiarity[participant_id],
1 - Familiarity[participant_id]) *
c(1 - Big5[participant_id],
1 - Big5[participant_id],
1 - Big5[participant_id]) *
c(1 - expressibility_score,
1 - expressibility_score,
1 - expressibility_score)
adjusted_prob <- adjusted_prob / sum(adjusted_prob)
n_attempts <- sample(1:3, 1, prob = adjusted_prob)
prev_answer_similarity <- NA # First attempt has no previous similarity
# Generate rows for each communicative attempt
for (attempt in 1:n_attempts) {
Eff <- base_effort # Start with base effort
# Modify effort for second and third attempts
if (attempt == 2) {
Eff <- Eff * b_comatt2
} else if (attempt == 3) {
Eff <- Eff * b_comatt3
}
# Adjust effort based on previous answer similarity
if (!is.na(prev_answer_similarity)) {
Eff <- Eff * (1 + (1 - prev_answer_similarity) * b_prevan)
}
# Store row
participant_data <- rbind(participant_data, data.frame(
Participant = participant_id,
Concept = concept_id,
Modality = modality,
Big5 = Big5[participant_id],
Familiarity = Familiarity[participant_id],
Expressibility = expressibility_score,
CommAtt = attempt,
Eff = Eff,
TrialNumber = trial_number,
PrevAn = prev_answer_similarity
))
# Update for next attempt
trial_number <- trial_number + 1
prev_answer_similarity <- runif(1, min = 0, max = 1) # Simulate next similarity
}
}
return(participant_data)
}
# Simulate data for all participants
final_data_synt <- do.call(rbind, lapply(participants, simulate_participant))
# Preview results
head(final_data_synt)## Participant Concept Modality Big5 Familiarity Expressibility CommAtt
## 1 1 1 combined 1.919899 1.843906 0.5761723 1
## 2 1 1 combined 1.919899 1.843906 0.5761723 2
## 3 1 1 combined 1.919899 1.843906 0.5761723 3
## 4 1 2 gesture 1.919899 1.843906 0.3493787 1
## 5 1 2 gesture 1.919899 1.843906 0.3493787 2
## 6 1 2 gesture 1.919899 1.843906 0.3493787 3
## Eff TrialNumber PrevAn
## 1 4.343849 1 NA
## 2 6.566250 2 0.98450637
## 3 2.367010 3 0.82035722
## 4 5.497853 4 NA
## 5 10.367605 5 0.48565959
## 6 3.998448 6 0.09090184So now we have synthetic data where we exactly defined what is the relationship between certain variables
- CommAtt -> Eff
The effort for second attempt is multiplied by 1.50 (Beta = 1.50) The effort for third attempts by 0.7 (ie decreases, Beta = 0.5)
- Fam -> Eff & CommAtt
Beta = 1.10 for Eff Negative effect on CommAtt
- Big5 -> Eff & CommAtt
Beta = 1.15 for Eff Negative effect on CommAtt
- Expr -> Eff & CommAtt
Beta = 1.20 for Eff Negative effect on CommAtt
(For simplicity reasons, we do not create an affect of trial number, as it does not represent a crucial variable)
Note that the synthetic data do not copy the structure of the data and the experimental design perfectly, but it does provide a reliable ground to build the models and test their performance.
Exploring structure
nrow(final_data_synt) # this is the number of datapoints## [1] 5076hist(final_data_synt$Eff)
final_data_synt |>
janitor::tabyl(Participant, Concept) # the number of repetition per concept by participant## Participant 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21
## 1 3 3 3 2 1 3 1 3 3 1 1 1 1 2 3 1 3 2 2 2 3
## 2 3 2 1 2 2 3 3 1 3 2 3 1 1 3 3 2 2 1 1 2 1
## 3 1 2 2 3 2 3 2 3 2 1 3 2 3 1 1 2 1 3 3 2 1
## 4 3 3 1 3 3 3 1 1 1 1 1 3 2 2 2 1 3 1 1 1 3
## 5 2 2 1 3 3 2 3 3 2 1 1 1 2 1 2 2 3 1 3 1 3
## 6 3 2 2 3 1 1 3 1 1 2 1 2 1 3 3 3 2 3 1 1 1
## 7 1 1 1 3 2 2 3 1 3 3 2 1 3 3 3 3 2 1 3 2 1
## 8 2 2 2 2 1 3 1 2 2 1 1 3 1 3 3 2 1 3 2 1 3
## 9 2 1 3 2 2 3 3 1 1 1 2 3 3 3 2 3 1 3 1 2 2
## 10 2 3 2 1 3 3 1 1 1 3 1 1 3 1 3 3 3 1 1 3 1
## 11 3 3 3 1 2 1 3 1 1 2 1 3 2 2 2 1 3 3 3 3 1
## 12 2 2 3 2 3 3 2 1 2 3 1 2 2 3 2 2 2 2 1 3 1
## 13 2 2 3 3 2 3 3 3 2 3 3 1 3 2 3 1 2 1 1 2 3
## 14 2 1 3 2 2 3 2 2 3 3 2 3 1 2 2 3 3 3 3 3 1
## 15 2 3 1 2 1 3 1 1 3 1 3 3 1 2 1 1 3 2 1 2 3
## 16 1 3 2 2 1 3 1 3 3 1 2 2 1 1 2 2 2 1 1 3 3
## 17 1 1 2 1 3 3 3 1 3 3 3 3 3 2 3 3 1 1 3 3 1
## 18 1 1 1 1 2 1 1 1 3 1 2 1 2 3 2 2 2 2 1 1 2
## 19 2 1 2 3 2 1 1 1 1 3 1 1 3 3 3 3 3 3 1 3 2
## 20 2 2 1 3 1 1 3 3 2 1 3 3 2 2 3 2 1 3 2 2 2
## 21 2 2 3 1 1 1 1 3 2 1 2 3 2 3 1 3 3 2 3 1 3
## 22 3 1 1 1 3 3 2 1 2 3 1 2 1 2 1 1 2 3 3 1 1
## 23 2 2 1 1 1 2 2 1 3 2 2 2 3 3 1 3 2 3 2 1 2
## 24 3 2 1 3 3 1 3 1 2 3 3 2 3 2 1 2 1 2 3 3 1
## 25 2 2 1 1 3 1 3 2 1 3 2 3 3 2 3 1 3 1 1 3 3
## 26 2 2 2 2 1 1 2 1 2 2 1 1 2 3 3 3 1 2 1 2 3
## 27 3 2 3 1 3 1 2 1 1 3 2 3 1 3 2 1 1 2 1 3 2
## 28 2 2 1 3 2 1 1 1 3 2 1 2 3 3 3 1 3 1 2 1 1
## 29 2 3 2 2 2 3 3 2 1 2 1 2 3 2 2 3 3 2 3 2 2
## 30 1 3 3 1 2 3 3 2 1 2 3 1 3 1 1 1 2 2 3 3 3
## 31 1 3 3 3 3 3 2 3 1 3 1 2 2 2 1 2 3 3 2 2 3
## 32 3 2 2 1 3 2 2 1 3 1 1 2 2 3 2 3 1 1 3 2 2
## 33 2 1 2 3 1 1 1 1 2 2 3 2 2 1 1 3 1 3 2 2 1
## 34 1 1 3 1 1 2 2 2 1 3 3 1 1 3 2 3 2 2 3 1 3
## 35 1 3 3 3 1 1 2 3 1 2 1 2 3 2 1 1 3 1 1 1 3
## 36 2 2 3 2 3 3 3 1 1 1 3 1 3 1 3 3 3 3 1 3 1
## 37 1 1 1 3 1 2 2 1 3 3 3 2 3 1 2 2 3 3 2 1 1
## 38 3 1 3 3 1 3 3 1 3 1 2 1 3 1 3 1 3 3 3 2 1
## 39 3 1 1 3 2 1 2 2 1 2 1 2 1 1 2 3 1 1 1 2 1
## 40 1 3 1 1 3 2 1 2 1 2 2 3 3 3 3 2 1 1 2 2 1
## 41 3 1 3 1 1 2 2 1 3 2 2 1 3 2 1 3 3 3 2 1 3
## 42 1 1 2 3 2 3 3 2 3 2 3 3 2 2 1 3 2 2 1 2 1
## 43 1 1 2 3 1 2 1 1 3 2 2 3 1 1 2 3 1 1 2 1 2
## 44 2 1 1 1 2 3 1 1 1 3 3 2 3 2 3 3 3 3 3 2 2
## 45 2 2 3 3 3 1 3 2 3 3 2 1 3 1 1 3 1 1 2 1 2
## 46 1 1 3 1 1 2 1 2 1 1 1 1 1 1 1 1 3 3 2 3 2
## 47 2 1 3 2 1 2 1 1 1 3 1 3 2 1 1 3 3 2 1 2 3
## 48 2 3 2 3 1 2 3 3 1 1 1 1 2 2 3 2 2 2 3 3 3
## 49 2 1 1 3 3 3 2 2 2 3 1 1 2 3 1 1 2 1 3 2 3
## 50 2 1 2 2 1 3 3 3 3 3 2 3 2 3 2 2 1 3 3 2 3
## 51 1 2 2 3 2 2 2 3 3 1 3 3 2 3 3 1 1 3 3 3 3
## 52 3 1 2 2 1 3 1 1 1 2 2 3 3 3 1 3 2 3 2 1 2
## 53 1 3 2 3 3 1 2 2 2 2 1 3 2 1 1 1 1 3 2 1 3
## 54 3 3 3 1 3 3 1 1 3 1 2 2 1 3 3 1 2 3 3 3 1
## 55 1 1 2 2 1 3 2 1 1 1 3 2 2 3 2 2 1 1 2 2 3
## 56 3 2 2 3 2 3 3 1 2 3 2 2 2 3 3 1 2 1 1 1 1
## 57 3 3 3 1 3 2 2 1 2 2 2 2 3 3 3 3 3 3 3 1 1
## 58 3 2 1 3 2 2 3 1 3 1 2 2 3 3 1 3 1 3 2 3 1
## 59 3 1 2 3 2 1 3 1 3 2 2 2 2 3 1 1 1 1 1 1 2
## 60 2 1 1 2 2 3 2 2 3 1 3 1 2 1 3 1 3 3 3 3 2
## 61 2 3 3 3 2 2 2 1 3 2 1 3 3 3 3 1 3 2 2 2 1
## 62 3 1 2 3 1 3 2 3 1 3 2 1 3 1 1 1 1 2 3 1 1
## 63 2 2 1 3 3 1 3 2 1 2 2 3 2 1 1 1 3 2 3 2 2
## 64 2 2 3 2 3 3 2 1 1 3 1 1 3 3 2 3 3 2 2 1 3
## 65 1 1 1 2 2 1 1 2 2 1 3 3 1 2 2 3 2 2 2 3 3
## 66 1 3 1 1 3 2 1 2 2 3 1 3 1 3 1 3 3 3 2 2 2
## 67 3 2 1 1 2 3 1 2 1 2 2 3 1 3 3 3 2 1 3 3 2
## 68 1 1 3 2 3 2 2 3 2 2 2 2 3 2 1 2 1 3 3 1 2
## 69 2 2 1 2 2 1 3 1 2 3 1 3 2 3 3 2 3 2 2 1 2
## 70 2 1 3 1 1 1 1 2 3 2 3 3 3 2 2 1 2 3 2 3 2
## 71 2 2 2 3 1 1 2 3 3 1 1 2 3 1 3 3 3 2 1 2 1
## 72 3 1 3 2 1 1 2 1 3 3 1 2 1 3 1 2 2 2 1 3 1
## 73 2 3 1 1 1 3 2 2 3 2 3 3 3 2 1 2 1 2 3 3 1
## 74 3 3 1 2 2 3 2 1 1 1 1 2 1 2 1 2 3 3 2 3 3
## 75 1 1 2 3 3 3 3 2 1 2 3 3 3 3 3 1 3 3 1 1 1
## 76 2 1 3 3 3 3 3 1 3 2 1 3 2 3 2 2 3 2 1 2 2
## 77 1 3 2 1 2 2 1 3 3 1 2 2 2 2 2 2 2 1 2 1 1
## 78 3 1 3 3 2 3 1 2 3 2 2 2 2 3 2 1 1 1 2 2 2
## 79 3 2 1 2 1 3 1 2 3 3 2 1 2 2 1 3 3 3 2 2 2
## 80 2 2 1 1 1 1 2 3 2 3 3 3 2 1 1 1 2 3 1 3 2
## 81 2 2 2 1 3 2 3 1 1 3 2 2 2 3 2 3 3 1 1 3 3
## 82 1 2 2 2 3 3 3 3 1 3 1 2 3 2 2 3 3 1 2 3 3
## 83 3 3 2 2 2 2 1 1 3 2 3 1 1 1 3 1 2 1 1 3 1
## 84 3 1 1 3 1 1 1 1 3 1 1 2 3 2 1 1 1 3 3 1 3
## 85 2 3 1 3 1 1 1 1 3 3 3 3 1 3 3 2 2 2 2 1 2
## 86 3 2 3 2 1 3 3 2 3 3 3 3 2 1 1 3 1 1 1 1 1
## 87 2 2 3 2 1 2 2 1 3 3 3 3 1 1 1 3 1 3 3 3 2
## 88 2 1 1 2 1 1 1 2 2 3 2 1 1 2 1 1 3 1 3 3 1
## 89 1 2 3 3 3 2 3 2 2 3 1 1 3 3 3 2 2 3 1 1 3
## 90 1 3 2 2 3 3 1 1 3 2 2 3 2 2 2 2 1 1 1 1 3
## 91 2 3 1 3 1 2 2 1 2 2 2 3 1 3 1 3 1 3 2 1 1
## 92 2 2 2 3 2 3 2 2 2 1 2 1 2 1 2 2 2 1 1 1 3
## 93 2 1 1 2 2 3 1 3 1 1 3 3 1 3 1 2 2 1 3 3 3
## 94 2 2 1 2 3 2 3 3 1 2 1 3 2 2 1 2 1 3 1 3 3
## 95 2 1 1 1 3 2 3 3 3 3 1 1 2 3 3 3 3 2 3 3 3
## 96 3 2 1 1 2 1 1 2 3 1 1 2 1 1 2 2 3 2 2 2 2
## 97 1 2 1 2 3 1 1 1 1 2 3 3 3 3 2 3 2 3 1 2 2
## 98 3 1 1 1 2 2 2 1 1 2 2 2 1 1 1 1 2 1 2 1 2
## 99 3 1 3 1 3 1 3 3 3 3 2 3 2 3 2 3 3 2 2 2 2
## 100 3 3 3 1 3 3 3 3 2 2 2 2 3 2 2 2 2 1 2 3 1
## 101 1 1 2 2 1 2 3 1 1 2 1 3 2 1 1 3 2 2 3 1 2
## 102 3 3 3 1 2 2 2 3 1 2 3 3 2 3 3 2 2 1 1 2 3
## 103 2 2 1 2 1 2 1 3 2 2 3 2 3 1 1 1 1 1 2 3 1
## 104 1 1 1 3 2 2 2 2 2 3 3 2 3 1 2 2 1 2 3 3 3
## 105 3 1 2 2 1 2 3 2 2 2 2 2 2 2 1 1 2 2 1 1 1
## 106 3 1 2 1 1 2 1 3 1 3 2 1 1 3 1 3 2 1 1 1 2
## 107 3 3 3 2 2 2 3 1 2 1 3 1 1 3 2 2 1 2 1 3 3
## 108 3 1 3 3 2 1 2 3 1 2 2 3 1 3 1 3 3 2 3 2 2
## 109 2 3 3 1 2 3 3 2 3 3 1 2 2 2 2 1 2 2 1 3 2
## 110 2 2 2 2 1 3 1 3 2 3 3 2 3 3 1 1 1 2 2 3 3
## 111 1 1 2 1 3 2 1 1 2 1 2 1 1 2 2 2 1 2 3 1 3
## 112 3 2 3 3 2 1 3 1 1 1 3 2 1 3 2 1 3 2 1 3 3
## 113 3 3 1 2 3 3 1 2 1 1 3 3 2 2 3 3 3 2 3 1 3
## 114 2 1 3 2 1 1 2 1 1 3 3 1 1 1 3 1 1 3 3 3 2
## 115 3 1 3 2 3 1 3 3 3 3 3 1 1 2 3 3 2 1 3 2 2
## 116 1 1 1 1 2 2 3 2 3 1 1 3 2 2 3 2 2 3 2 1 1
## 117 2 3 2 3 2 1 1 3 1 1 2 2 3 2 3 3 3 1 2 1 1
## 118 3 2 1 2 2 3 1 1 1 3 2 2 2 1 3 1 3 3 1 2 3
## 119 1 3 1 2 2 2 1 2 1 2 2 3 3 1 2 1 3 3 2 3 1
## 120 1 1 3 1 1 3 2 1 3 3 2 2 1 3 3 2 1 2 3 1 1Plot: H1
# Create a boxplot comparing Effort across different Communicative Attempts
ggplot(final_data_synt, aes(x = as.factor(CommAtt), y = Eff)) +
geom_boxplot(aes(fill = as.factor(CommAtt))) +
labs(title = "Comparison of Effort Across Communicative Attempts",
x = "Communicative Attempts",
y = "Effort",
fill = "CommAtt") +
theme_minimal() +
theme(legend.position = "none") 
The visual already hints on the effects we expect to exist in the
data.
Plot: H2
# Filter out CommAtt == 1
filtered_data <- final_data_synt[final_data_synt$CommAtt != 1, ]
# Scatter plot with regression line
library(ggplot2)
ggplot(filtered_data, aes(x = PrevAn, y = Eff)) +
geom_point(alpha = 0.6, color = "blue") + # Scatter points
geom_smooth(method = "lm", color = "red", se = FALSE) + # Regression line
labs(x = "Previous Answer Similarity (PrevAn)",
y = "Effort (Eff)",
title = "Relationship between Effort and Previous Answer Similarity") +
theme_minimal()
H1: Stating causal model
Before doing that, we create a causal diagram, so-called directed acyclic graphs (DAG). This will help us clarify the causal model and introduce the predictors we aim to model.
Our two hypothesis go as follows:
H1: Correction recruits more physical effort than the baseline performance.
H2: A higher degree of misunderstanding will require a performer to engage in more effortful correction.
The relationship of interest is the causal effect of communicative attempt on effort. Our assumptions include:
Personality traits (measured with Big5) will influence effort (e.g., people are more willing to put more effort if they are open-minded) and also communicative attempt (e.g., more extroverted people are better in this game, therefore they need less attempts)
Familiarity with guessing partner influences effort (ditto) as well as communicative attempt (e.g., friends are better in this game than strangers)
Similarly, participant will also directly influence effort and commAtt, because personality traits might not be capturing all variation
Expressibility of concepts is going to influence effort (e.g., more expressible concepts allow more enhancement - but could be also other direction) as well as CommAtt (e.g., more expressible concepts are more readily guessed and don’t need more attempts)
Similarly, concept will also directly influence effort and commAtt, because expressibility might not be capturing all variation
Modality (uni vs. multi) will influence the value of effort. We assume that in unimodal condition, the feature does not need to account for synergic relations with the other modality, and may carry the whole signal. In multimodal condition, however, there may be synergic relations that moderate the full expressiveness of this feature
Trial number (characterising how far one is in the experiment) could be hinting on learning processes through out the experiment, or - the other direction - on increasing fatigue
Previous answer (PrevAn) will affect the effort (more similar answer will require less effortful correction) as well as communicative attempt (correct answers do not require further corrections)
Expressibility of concepts will influence the answer (more expressible are easier to guess)
Trial number will also affect previous answer if there is learning involved, or conversely, increasing fatigue
These are the implied conditional independencies
## Big5 _||_ Conc
## Big5 _||_ Expr
## Big5 _||_ Fam | Pcn
## Big5 _||_ Mod
## Big5 _||_ TrNm
## Conc _||_ Fam
## Conc _||_ Mod
## Conc _||_ Pcn
## Conc _||_ TrNm
## Expr _||_ Fam
## Expr _||_ Mod
## Expr _||_ Pcn
## Expr _||_ TrNm
## Fam _||_ Mod
## Fam _||_ TrNm
## Mod _||_ Pcn
## Mod _||_ TrNm
## Pcn _||_ TrNmThis is the adjustment set that needs to be included in the model to make sure we are not confounding
## { Big5, Conc, Expr, Fam, Mod, Pcn, TrNum }H1: Modelling
In this part, we are going to test the following hypothesis:
H1: Correction recruits more physical effort than the baseline performance.
Contrast coding
To be able to interpret the data, we will first need to do some contrast coding
- Convert columns to factors
final_data$CommAtt <- as.factor(final_data$CommAtt)
final_data$Modality <- as.factor(final_data$Modality)
final_data$Participant <- as.factor(final_data$Participant)
final_data$Concept <- as.factor(final_data$Concept)
final_data$TrialNumber <- as.numeric(final_data$TrialNumber) # Ensure TrialNumber is numeric- Contrast-coding of categorical variables
contrasts(final_data$CommAtt) <- MASS::contr.sdif(3)
contrasts(final_data$Modality) <- contr.sum(3)/2This is how CommAtt is contrast-coded now
2-1 3-21 -0.6666667 -0.3333333 2 0.3333333 -0.3333333 3 0.3333333 0.6666667
This is how modality is cc-ed
[,1] [,2]combined 0.5 0.0 gesture 0.0 0.5 vocal -0.5 -0.5
- Center continuous variables
final_data$TrialNumber_c <- final_data$TrialNumber - median(range(final_data$TrialNumber))
final_data$Familiarity <- final_data$Familiarity - median(range(final_data$Familiarity))
final_data$Big5 <- final_data$Big5 - median(range(final_data$Big5))
# For now, we will just center Familiarity and Big5 because we created them synthetically. But the real data have these two variables already z-scored- Z-score expressibility within a modality
final_data <-
final_data |>
group_by(Modality) |>
mutate(Expressibility_z = (Expressibility - mean(Expressibility))/ sd(final_data$Expressibility, na.rm = T)) |>
ungroup()Model 1 - Simple reproduction of DAG
First we want to do a simple model that, however, reproduce our DAG.
Our main predictor is communicative attempt (CommAtt). To account for confounders - variables affecting both the predictor and the dependent variable - we need to adjust for them in the model by including them as covariates to isolate the effect of the predictor. Based on our DAG, confounders include:
- familiarity
- big5
- expressibility
- trial number
- modality
- concept
- participant
We include 1.-5. as fixed factors. For 6.-7., we include varying intercepts as we expect that each participant and concept may have they own baseline level of effort and thus allow for individual variation. Partial pooling is also beneficial in that extreme values (or categories will fewer datapoints) will be shrunk toward the overal average.
We will now not include PrevAn (previous answer) variable because we will need to do some further data-wrangling when building model for H2. That is mainly because PrevAn has some NA values, concretely for first attempts. Models would automatically exclude NA data, and we would therefore loose all effort data for attempt 1. For H1, however, we want to keep it there.
This is OUR FIRST model without setting any priors, leaving them to default values
h1.m1 <- brm(Eff ~ 1 + CommAtt + Familiarity + Big5 + Expressibility_z + TrialNumber_c + Modality + (1 | Participant) + (1 | Concept),
data = final_data,
iter = 4000,
cores = 4)
# Add criterions for later diagnostics
h1.m1 <- add_criterion(h1.m1, criterion = c("loo", "waic"))
# Calculate also variance explained (R^2)
h1.m1_R2 <- bayes_R2(h1.m1)
# Save both as objects
saveRDS(h1.m1, here("09_Analysis_Modeling", "models", "h1.m1.rds"))
saveRDS(h1.m1_R2, here("09_Analysis_Modeling", "models", "h1.m1_R2.rds"))
beep(5)h1.m1 <- readRDS(here("09_Analysis_Modeling", "models", "h1.m1.rds"))
h1.m1_R2 <- readRDS(here("09_Analysis_Modeling", "models", "h1.m1_R2.rds"))
# Summary
summary(h1.m1)## Family: gaussian
## Links: mu = identity; sigma = identity
## Formula: Eff ~ 1 + CommAtt + Familiarity + Big5 + Expressibility_z + TrialNumber_c + Modality + (1 | Participant) + (1 | Concept)
## Data: final_data (Number of observations: 5076)
## Draws: 4 chains, each with iter = 4000; warmup = 2000; thin = 1;
## total post-warmup draws = 8000
##
## Multilevel Hyperparameters:
## ~Concept (Number of levels: 21)
## Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
## sd(Intercept) 0.04 0.02 0.00 0.09 1.00 2122 1863
##
## ~Participant (Number of levels: 120)
## Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
## sd(Intercept) 0.03 0.02 0.00 0.07 1.00 2320 3475
##
## Regression Coefficients:
## Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
## Intercept 4.02 0.02 3.98 4.06 1.00 6694 5622
## CommAtt2M1 3.08 0.03 3.03 3.14 1.00 8963 6147
## CommAtt3M2 -4.39 0.04 -4.46 -4.32 1.00 8289 5880
## Familiarity 1.26 0.02 1.21 1.30 1.00 8380 5226
## Big5 1.29 0.02 1.25 1.34 1.00 9979 6018
## Expressibility_z 0.46 0.02 0.43 0.49 1.00 9888 5989
## TrialNumber_c -0.00 0.00 -0.00 0.00 1.00 7595 5984
## Modality1 -1.30 0.04 -1.38 -1.23 1.00 8559 5997
## Modality2 0.64 0.04 0.57 0.71 1.00 7347 5787
##
## Further Distributional Parameters:
## Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
## sigma 0.91 0.01 0.89 0.93 1.00 11791 5667
##
## Draws were sampled using sampling(NUTS). For each parameter, Bulk_ESS
## and Tail_ESS are effective sample size measures, and Rhat is the potential
## scale reduction factor on split chains (at convergence, Rhat = 1).# Note on modality: Modality1 represents the baseline difference between combined and the other modalities, and Modality2 represents the difference between gesture and vocal compared to combined.Transforming coefficients
To be able to link these estimates back to the simulation coefficients, let’s create a function
transform_attempt <- function(intercept, CommAtt2M1, CommAtt3M2) {
# Effort for the first attempt (base effort)
Eff_attempt_1 <- intercept
# Effort for the second attempt
Eff_attempt_2 <- Eff_attempt_1 + CommAtt2M1
# Effort for the third attempt
Eff_attempt_3 <- Eff_attempt_1 + CommAtt2M1 + CommAtt3M2
# Calculate ratios
b_attempt2 <- Eff_attempt_2 / Eff_attempt_1
b_attempt3 <- Eff_attempt_3 / Eff_attempt_2
return(data.frame(b_attempt2, b_attempt3))
}
h1m1_coeff <- transform_attempt(4.02, 3.08, -4.39)
print(h1m1_coeff)## b_attempt2 b_attempt3
## 1 1.766169 0.3816901# For centered familiarity
fam = (4.02 + 1.26) / 4.02
print(fam)## [1] 1.313433# For centered BIF
bif = (4.02 + 1.29) / 4.02
print(bif)## [1] 1.320896# Expressibility is z-scored so we will not get to the coefficient in the same way but we can check the conditinal effects for checking whether it looks goodLet’s also check the visuals
plot(h1.m1)
# all caterpillars look nice
plot(conditional_effects(h1.m1), points = TRUE)
# the effects all go in good direction
pp_check(h1.m1, type = "dens_overlay")
# Looks good but not amazing - mostly because the posteriors seem to not know effort cannot be negative
pp_check(h1.m1, type = "error_scatter_avg")
# half blobby, half correlated, so still some room for improvement
# positive correlation means that errors increase with predicted values. So the model does perform well for some range, but becomes less reliable with increase in the predicted values
# also the blob is centered around 0 which is good
# it could be we are ignoring some interaction terms or non-linearity (which we know we kind of do). Transformation could also help (e.g., log). Of course, we are also still not specifying any priors so let's not yet make it a disaster
h1.m1_R2## Estimate Est.Error Q2.5 Q97.5
## R2 0.8381496 0.001662319 0.8348364 0.8412852# explained variance around 83%Overall, we see good directions of all predictors, mostly also in accordance with the expected coefficients. Of course, the synthetic data is quite complex so there might be other dependencies that moderate the causal relationships and that is why we do not see exactly the numbers we use to create the data.
Let’s have another model for comparison.
The fixed predictors make sense like this, as it will allow us to assess the effect of each on the effort, despite it not being our main research question. But we can assume that participants and concept have not only different baselines of effort (varying intercept). The effect of CommAtt on effort might vary across them too, hence we can try to add varying slopes for them and see whether the diagnostics improves. We will also add TrialNumber as a varying intercept, because we expect variation between earlier and later performances (because of learning, or opposite, fatigue) and we do not really need a single coefficient for this predictor anyway.
Model 2 - Varying slopes and intercepts
h1.m2 <- brm(Eff ~ 1 + CommAtt + Modality + Big5 + Familiarity + Expressibility_z + (1 + CommAtt | Participant) + (1 + CommAtt | Concept) + (1 | TrialNumber_c),
data = final_data,
iter = 4000,
cores = 4)
# Add criterions for later diagnostics
h1.m2 <- add_criterion(h1.m2, criterion = c("loo", "waic"))
# Calculate also variance explained (R^2)
h1.m2_R2 <- bayes_R2(h1.m2)
# Save both as objects
saveRDS(h1.m2, here("09_Analysis_Modeling", "models", "h1.m2.rds"))
saveRDS(h1.m2_R2, here("09_Analysis_Modeling", "models", "h1.m2_R2.rds"))
beep(5)h1.m2 <- readRDS(here("09_Analysis_Modeling", "models", "h1.m2.rds"))
h1.m2_R2 <- readRDS(here("09_Analysis_Modeling", "models", "h1.m2_R2.rds"))
# Summary
summary(h1.m2)## Family: gaussian
## Links: mu = identity; sigma = identity
## Formula: Eff ~ 1 + CommAtt + Modality + Big5 + Familiarity + Expressibility_z + (1 + CommAtt | Participant) + (1 + CommAtt | Concept) + (1 | TrialNumber_c)
## Data: final_data (Number of observations: 5076)
## Draws: 4 chains, each with iter = 4000; warmup = 2000; thin = 1;
## total post-warmup draws = 8000
##
## Multilevel Hyperparameters:
## ~Concept (Number of levels: 21)
## Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS
## sd(Intercept) 0.06 0.02 0.02 0.10 1.00 2791
## sd(CommAtt2M1) 0.17 0.04 0.11 0.26 1.00 3239
## sd(CommAtt3M2) 0.25 0.06 0.16 0.37 1.00 3386
## cor(Intercept,CommAtt2M1) 0.52 0.28 -0.14 0.93 1.00 1312
## cor(Intercept,CommAtt3M2) -0.35 0.32 -0.90 0.32 1.00 1261
## cor(CommAtt2M1,CommAtt3M2) -0.93 0.07 -1.00 -0.75 1.00 5691
## Tail_ESS
## sd(Intercept) 2048
## sd(CommAtt2M1) 4832
## sd(CommAtt3M2) 5133
## cor(Intercept,CommAtt2M1) 1617
## cor(Intercept,CommAtt3M2) 1450
## cor(CommAtt2M1,CommAtt3M2) 6716
##
## ~Participant (Number of levels: 120)
## Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS
## sd(Intercept) 0.15 0.04 0.08 0.22 1.00 901
## sd(CommAtt2M1) 0.78 0.06 0.68 0.90 1.00 2066
## sd(CommAtt3M2) 1.10 0.08 0.95 1.27 1.00 2076
## cor(Intercept,CommAtt2M1) 0.96 0.03 0.87 1.00 1.01 352
## cor(Intercept,CommAtt3M2) -0.96 0.04 -1.00 -0.85 1.01 353
## cor(CommAtt2M1,CommAtt3M2) -1.00 0.00 -1.00 -0.99 1.00 7139
## Tail_ESS
## sd(Intercept) 2511
## sd(CommAtt2M1) 3409
## sd(CommAtt3M2) 3623
## cor(Intercept,CommAtt2M1) 596
## cor(Intercept,CommAtt3M2) 690
## cor(CommAtt2M1,CommAtt3M2) 6841
##
## ~TrialNumber_c (Number of levels: 50)
## Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
## sd(Intercept) 0.02 0.02 0.00 0.06 1.00 3119 3942
##
## Regression Coefficients:
## Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
## Intercept 4.03 0.02 3.99 4.08 1.00 2093 3494
## CommAtt2M1 3.07 0.08 2.90 3.23 1.00 1384 2570
## CommAtt3M2 -4.36 0.12 -4.59 -4.13 1.00 1411 2477
## Modality1 -1.29 0.03 -1.35 -1.23 1.00 8986 6966
## Modality2 0.64 0.03 0.58 0.71 1.00 9767 6458
## Big5 1.08 0.05 0.99 1.16 1.01 854 2264
## Familiarity 1.04 0.05 0.94 1.12 1.00 957 2173
## Expressibility_z 0.44 0.02 0.41 0.47 1.00 5445 5281
##
## Further Distributional Parameters:
## Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
## sigma 0.80 0.01 0.78 0.81 1.00 11536 6199
##
## Draws were sampled using sampling(NUTS). For each parameter, Bulk_ESS
## and Tail_ESS are effective sample size measures, and Rhat is the potential
## scale reduction factor on split chains (at convergence, Rhat = 1).# Coefficients remain mostly unchanged but there were some divergent transitionsplot(h1.m2)
# some some of the caterpillars are not that pretty anymore
plot(conditional_effects(h1.m2), points = TRUE)
# the effects all go in good direction
pp_check(h1.m2, type = "dens_overlay")
# Looks good but not amazing - mostly because the posteriors seem to not know effort cannot be negative
pp_check(h1.m2, type = "error_scatter_avg")
# This looks a bit more blobby than before but still lots of errors for higher values
h1.m2_R2## Estimate Est.Error Q2.5 Q97.5
## R2 0.8760573 0.00125936 0.8734939 0.8784314# explained variance around 87%One of the reasons why we might be getting the divergent transition is because of the correlation between slopes and intercepts in the previous model. Let’s now get rid of it
Model 3 - No correlation coefficient
h1.m3 <- brm(Eff ~ 1 + CommAtt + Modality + Big5 + Familiarity + Expressibility_z + (1 + CommAtt || Participant) + (1 + CommAtt || Concept) + (1 || TrialNumber_c),
data = final_data,
iter = 4000,
cores = 4)
# Add criterions for later diagnostics
h1.m3 <- add_criterion(h1.m3, criterion = c("loo", "waic"))
# Calculate also variance explained (R^2)
h1.m3_R2 <- bayes_R2(h1.m3)
# Save both as objects
saveRDS(h1.m3, here("09_Analysis_Modeling", "models", "h1.m3.rds"))
saveRDS(h1.m3_R2, here("09_Analysis_Modeling", "models", "h1.m3_R2.rds"))
beep(5)h1.m3 <- readRDS(here("09_Analysis_Modeling", "models", "h1.m3.rds"))
h1.m3_R2 <- readRDS(here("09_Analysis_Modeling", "models", "h1.m3_R2.rds"))
# Summary
summary(h1.m3)## Family: gaussian
## Links: mu = identity; sigma = identity
## Formula: Eff ~ 1 + CommAtt + Modality + Big5 + Familiarity + Expressibility_z + (1 + CommAtt || Participant) + (1 + CommAtt || Concept) + (1 || TrialNumber_c)
## Data: final_data (Number of observations: 5076)
## Draws: 4 chains, each with iter = 4000; warmup = 2000; thin = 1;
## total post-warmup draws = 8000
##
## Multilevel Hyperparameters:
## ~Concept (Number of levels: 21)
## Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
## sd(Intercept) 0.05 0.02 0.01 0.10 1.00 2648 2360
## sd(CommAtt2M1) 0.09 0.05 0.01 0.19 1.01 1159 1540
## sd(CommAtt3M2) 0.17 0.06 0.04 0.30 1.00 1422 1517
##
## ~Participant (Number of levels: 120)
## Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
## sd(Intercept) 0.03 0.02 0.00 0.08 1.00 2442 4123
## sd(CommAtt2M1) 0.71 0.05 0.61 0.82 1.00 2733 4516
## sd(CommAtt3M2) 1.05 0.08 0.91 1.21 1.00 2231 3994
##
## ~TrialNumber_c (Number of levels: 50)
## Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
## sd(Intercept) 0.02 0.02 0.00 0.06 1.00 3301 3712
##
## Regression Coefficients:
## Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
## Intercept 4.03 0.02 3.99 4.06 1.00 9303 6039
## CommAtt2M1 3.07 0.07 2.93 3.22 1.00 2152 3434
## CommAtt3M2 -4.36 0.11 -4.58 -4.15 1.00 2100 3982
## Modality1 -1.30 0.03 -1.36 -1.24 1.00 11936 6326
## Modality2 0.65 0.03 0.58 0.71 1.00 12551 6401
## Big5 1.23 0.02 1.19 1.28 1.00 15911 6657
## Familiarity 1.20 0.02 1.15 1.24 1.00 14065 6574
## Expressibility_z 0.45 0.01 0.42 0.48 1.00 14842 5619
##
## Further Distributional Parameters:
## Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
## sigma 0.81 0.01 0.80 0.83 1.00 11557 5519
##
## Draws were sampled using sampling(NUTS). For each parameter, Bulk_ESS
## and Tail_ESS are effective sample size measures, and Rhat is the potential
## scale reduction factor on split chains (at convergence, Rhat = 1).# Coefficients remain mostly unchanged plot(h1.m3)
# no correlation fixed the performance
plot(conditional_effects(h1.m3), points = TRUE)
# the effects all go in good direction
pp_check(h1.m3, type = "dens_overlay")
# the problem with predicting negative values remains
pp_check(h1.m3, type = "error_scatter_avg")
# unchanged
h1.m3_R2## Estimate Est.Error Q2.5 Q97.5
## R2 0.8714492 0.001445024 0.8685275 0.874215# explained variance remains around 87%Model 3.1 - adding priors
priors_h1m3p <- c(
set_prior("normal(2.5, 0.5)", class = "Intercept", lb=0),
set_prior("normal(0,0.50)", class = "b", coef = "CommAtt2M1"),
set_prior("normal(0,0.50)", class = "b", coef = "CommAtt3M2"),
set_prior("normal(0,0.25)", class = "b", coef = "Modality1"),
set_prior("normal(0,0.25)", class = "b", coef = "Modality2"),
set_prior("normal(0,0.25)", class = "b", coef = "Big5"),
set_prior("normal(0,0.25)", class = "b", coef = "Familiarity"),
set_prior("normal(0,0.25)", class = "b", coef = "Expressibility_z"),
set_prior("normal(0.5,0.1)", class = "sd", group = "TrialNumber_c"),
set_prior("normal(0.5,0.1)", class = "sd", group = "Participant"),
set_prior("normal(0.5,0.1)", class = "sd", group = "Concept"),
set_prior("normal(1,0.1)", class = "sd"),
set_prior("normal(0.5,0.25)", class = "sigma")
)
h1.m3p <- brm(Eff ~ 1 + CommAtt + Modality + Big5 + Familiarity + Expressibility_z + (1 + CommAtt || Participant) + (1 + CommAtt || Concept) + (1 || TrialNumber_c),
data = final_data,
prior=priors_h1m3p,
family = gaussian,
iter = 4000,
cores = 4)
# Add criterions for later diagnostics
h1.m3p <- add_criterion(h1.m3p, criterion = c("loo", "waic"))
# Calculate also variance explained (R^2)
h1.m3p_R2 <- bayes_R2(h1.m3p)
# Save both as objects
saveRDS(h1.m3p, here("09_Analysis_Modeling", "models", "h1.m3p.rds"))
saveRDS(h1.m3p_R2, here("09_Analysis_Modeling", "models", "h1.m3p_R2.rds"))
beep(5)h1.m3p <- readRDS(here("09_Analysis_Modeling", "models", "h1.m3p.rds"))
h1.m3p_R2 <- readRDS(here("09_Analysis_Modeling", "models", "h1.m3p_R2.rds"))
# Summary
summary(h1.m3p)## Family: gaussian
## Links: mu = identity; sigma = identity
## Formula: Eff ~ 1 + CommAtt + Modality + Big5 + Familiarity + Expressibility_z + (1 + CommAtt || Participant) + (1 + CommAtt || Concept) + (1 || TrialNumber_c)
## Data: final_data (Number of observations: 5076)
## Draws: 4 chains, each with iter = 4000; warmup = 2000; thin = 1;
## total post-warmup draws = 8000
##
## Multilevel Hyperparameters:
## ~Concept (Number of levels: 21)
## Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
## sd(Intercept) 0.08 0.03 0.03 0.14 1.00 2732 4188
## sd(CommAtt2M1) 0.22 0.06 0.12 0.36 1.00 2853 4652
## sd(CommAtt3M2) 0.34 0.08 0.20 0.51 1.00 2590 5063
##
## ~Participant (Number of levels: 120)
## Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
## sd(Intercept) 0.06 0.02 0.01 0.10 1.00 2220 2200
## sd(CommAtt2M1) 0.66 0.04 0.58 0.75 1.00 4131 5851
## sd(CommAtt3M2) 0.89 0.05 0.80 0.98 1.00 4638 6332
##
## ~TrialNumber_c (Number of levels: 50)
## Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
## sd(Intercept) 0.04 0.02 0.00 0.08 1.00 2687 3691
##
## Regression Coefficients:
## Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
## Intercept 4.02 0.02 3.97 4.07 1.00 8469 6434
## CommAtt2M1 2.99 0.08 2.83 3.15 1.00 2795 4372
## CommAtt3M2 -4.12 0.12 -4.35 -3.87 1.00 2401 3679
## Modality1 -1.27 0.03 -1.34 -1.21 1.00 13767 6673
## Modality2 0.63 0.03 0.56 0.69 1.00 13131 6444
## Big5 1.23 0.02 1.18 1.28 1.00 12639 6019
## Familiarity 1.19 0.03 1.14 1.24 1.00 13068 6412
## Expressibility_z 0.45 0.02 0.42 0.48 1.00 15210 5985
##
## Further Distributional Parameters:
## Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
## sigma 0.81 0.01 0.80 0.83 1.00 14125 5732
##
## Draws were sampled using sampling(NUTS). For each parameter, Bulk_ESS
## and Tail_ESS are effective sample size measures, and Rhat is the potential
## scale reduction factor on split chains (at convergence, Rhat = 1).# Coefficients remain mostly unchanged plot(h1.m3p)
# all good
plot(conditional_effects(h1.m3p), points = TRUE)
# the effects all go in good direction
pp_check(h1.m3p, type = "dens_overlay")
# the problem with predicting negative values remains
pp_check(h1.m3p, type = "error_scatter_avg")
# unchanged
h1.m3p_R2## Estimate Est.Error Q2.5 Q97.5
## R2 0.8705705 0.001478901 0.8676255 0.873398# explained variance remains around 87%Model 4 - Restricting priors with exponential distribution
Let’s do one more test with the priors, restricting sd them with exponential distribution, i.e., negative values are not possible
priors_h1m4 <- c(
set_prior("normal(3, 0.3)", class = "Intercept", lb = 0),
set_prior("normal(0, 0.25)", class = "b", coef = "CommAtt2M1"),
set_prior("normal(0, 0.25)", class = "b", coef = "CommAtt3M2"),
set_prior("normal(0, 0.15)", class = "b", coef = "Modality1"),
set_prior("normal(0, 0.15)", class = "b", coef = "Modality2"),
set_prior("normal(0, 0.15)", class = "b", coef = "Big5"),
set_prior("normal(0, 0.15)", class = "b", coef = "Familiarity"),
set_prior("normal(0, 0.15)", class = "b", coef = "Expressibility_z"),
# Exponential priors for standard deviations
set_prior("exponential(3)", class = "sd", group = "TrialNumber_c"), # exp(3) has a mean of 1/3 and concentrates most density around small values
set_prior("exponential(3)", class = "sd", group = "Participant"),
set_prior("exponential(3)", class = "sd", group = "Concept"),
set_prior("exponential(1)", class = "sd"), # Generic sd prior
# Residual standard deviation - keep it narrow
set_prior("normal(0.5, 0.1)", class = "sigma")
)
h1.m4 <- brm(Eff ~ 1 + CommAtt + Modality + Big5 + Familiarity + Expressibility_z + (1 + CommAtt || Participant) + (1 + CommAtt || Concept) + (1 || TrialNumber_c),
data = final_data,
prior=priors_h1m4,
family = gaussian,
iter = 4000,
cores = 4)
# Add criterions for later diagnostics
h1.m4 <- add_criterion(h1.m4, criterion = c("loo", "waic"))
# Calculate also variance explained (R^2)
h1.m4_R2 <- bayes_R2(h1.m4)
# Save both as objects
saveRDS(h1.m4, here("09_Analysis_Modeling", "models", "h1.m4.rds"))
saveRDS(h1.m4_R2, here("09_Analysis_Modeling", "models", "h1.m4_R2.rds"))
beep(5)h1.m4 <- readRDS(here("09_Analysis_Modeling", "models", "h1.m4.rds"))
h1.m4_R2 <- readRDS(here("09_Analysis_Modeling", "models", "h1.m4_R2.rds"))
# Summary
summary(h1.m4)## Family: gaussian
## Links: mu = identity; sigma = identity
## Formula: Eff ~ 1 + CommAtt + Modality + Big5 + Familiarity + Expressibility_z + (1 + CommAtt || Participant) + (1 + CommAtt || Concept) + (1 || TrialNumber_c)
## Data: final_data (Number of observations: 5076)
## Draws: 4 chains, each with iter = 4000; warmup = 2000; thin = 1;
## total post-warmup draws = 8000
##
## Multilevel Hyperparameters:
## ~Concept (Number of levels: 21)
## Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
## sd(Intercept) 0.05 0.02 0.01 0.09 1.00 2376 2462
## sd(CommAtt2M1) 1.98 0.37 1.33 2.78 1.00 1362 2663
## sd(CommAtt3M2) 3.36 0.46 2.56 4.38 1.00 942 2031
##
## ~Participant (Number of levels: 120)
## Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
## sd(Intercept) 0.03 0.02 0.00 0.08 1.00 2217 4211
## sd(CommAtt2M1) 0.71 0.05 0.61 0.82 1.00 2746 4486
## sd(CommAtt3M2) 1.05 0.08 0.92 1.22 1.00 2446 4360
##
## ~TrialNumber_c (Number of levels: 50)
## Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
## sd(Intercept) 0.02 0.02 0.00 0.06 1.00 3383 3814
##
## Regression Coefficients:
## Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
## Intercept 4.02 0.02 3.99 4.06 1.00 9553 6582
## CommAtt2M1 0.88 0.30 0.31 1.47 1.00 4070 5104
## CommAtt3M2 -0.36 0.26 -0.87 0.14 1.00 10008 6008
## Modality1 -1.23 0.03 -1.29 -1.16 1.00 14208 6281
## Modality2 0.59 0.03 0.53 0.65 1.00 15010 6334
## Big5 1.21 0.02 1.16 1.25 1.00 13972 6575
## Familiarity 1.17 0.02 1.12 1.21 1.00 13947 5839
## Expressibility_z 0.44 0.01 0.41 0.47 1.00 17170 6516
##
## Further Distributional Parameters:
## Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
## sigma 0.81 0.01 0.79 0.83 1.00 15722 5841
##
## Draws were sampled using sampling(NUTS). For each parameter, Bulk_ESS
## and Tail_ESS are effective sample size measures, and Rhat is the potential
## scale reduction factor on split chains (at convergence, Rhat = 1).# the coefficients for commatt now seem much lowerh1m1_coeff <- transform_attempt(4.02, 0.88, -0.36)
print(h1m1_coeff)## b_attempt2 b_attempt3
## 1 1.218905 0.9265306# For centered familiarity
fam = (4.02 + 1.17) / 4.02
print(fam)## [1] 1.291045# For centered BIF
bif = (4.02 + 1.21) / 4.02
print(bif)## [1] 1.300995So it seems that we are misleading with the priors and this will not be the way
plot(h1.m4)
# all good
plot(conditional_effects(h1.m4), points = TRUE)
# the effect of main predictor is now very moderated
pp_check(h1.m4, type = "dens_overlay")
# the problem with predicting negative values remains
pp_check(h1.m4, type = "error_scatter_avg")
# unchanged
h1.m4_R2## Estimate Est.Error Q2.5 Q97.5
## R2 0.8701813 0.001472199 0.867155 0.872996# explained variance remains around 87%Model 5 - student family
So we still see negative values in the posterior simulations, so let’s try Student’s t-distribution which is more robut to outliers and can potentially reduce the likelihood of negative values (if we reduce degrees of freedom)
priors_h1m5 <- c(
set_prior("normal(2.5, 0.5)", class = "Intercept", lb=0),
set_prior("normal(0,0.50)", class = "b", coef = "CommAtt2M1"),
set_prior("normal(0,0.50)", class = "b", coef = "CommAtt3M2"),
set_prior("normal(0,0.25)", class = "b", coef = "Modality1"),
set_prior("normal(0,0.25)", class = "b", coef = "Modality2"),
set_prior("normal(0,0.25)", class = "b", coef = "Big5"),
set_prior("normal(0,0.25)", class = "b", coef = "Familiarity"),
set_prior("normal(0,0.25)", class = "b", coef = "Expressibility_z"),
set_prior("normal(0.5,0.05)", class = "sd", group = "TrialNumber_c"),
set_prior("normal(0.5,0.05)", class = "sd", group = "Participant"),
set_prior("normal(0.5,0.05)", class = "sd", group = "Concept"),
set_prior("normal(1,0.05)", class = "sd"),
set_prior("normal(0.5,0.1)", class = "sigma"),
set_prior("gamma(2, 0.1)", class = "nu") # Prior for degrees of freedom
)
h1.m5 <- brm(Eff ~ 1 + CommAtt + Modality + Big5 + Familiarity + Expressibility_z + (1 + CommAtt || Participant) + (1 + CommAtt || Concept) + (1 || TrialNumber_c),
data = final_data,
prior=priors_h1m5,
family = student,
iter = 4000,
cores = 4)
# Add criterions for later diagnostics
h1.m5 <- add_criterion(h1.m5, criterion = c("loo", "waic"))
# Calculate also variance explained (R^2)
h1.m5_R2 <- bayes_R2(h1.m5)
# Save both as objects
saveRDS(h1.m5, here("09_Analysis_Modeling", "models", "h1.m5.rds"))
saveRDS(h1.m5_R2, here("09_Analysis_Modeling", "models", "h1.m5_R2.rds"))
beep(5)h1.m5 <- readRDS(here("09_Analysis_Modeling", "models", "h1.m5.rds"))
h1.m5_R2 <- readRDS(here("09_Analysis_Modeling", "models", "h1.m5_R2.rds"))
# Summary
summary(h1.m5)## Family: student
## Links: mu = identity; sigma = identity; nu = identity
## Formula: Eff ~ 1 + CommAtt + Modality + Big5 + Familiarity + Expressibility_z + (1 + CommAtt || Participant) + (1 + CommAtt || Concept) + (1 || TrialNumber_c)
## Data: final_data (Number of observations: 5076)
## Draws: 4 chains, each with iter = 4000; warmup = 2000; thin = 1;
## total post-warmup draws = 8000
##
## Multilevel Hyperparameters:
## ~Concept (Number of levels: 21)
## Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
## sd(Intercept) 0.38 0.06 0.26 0.49 1.00 5974 5242
## sd(CommAtt2M1) 0.42 0.05 0.32 0.53 1.00 8447 6363
## sd(CommAtt3M2) 0.46 0.05 0.36 0.56 1.00 8880 6408
##
## ~Participant (Number of levels: 120)
## Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
## sd(Intercept) 0.13 0.02 0.10 0.17 1.00 3617 5313
## sd(CommAtt2M1) 0.57 0.03 0.51 0.63 1.00 6105 6010
## sd(CommAtt3M2) 0.72 0.03 0.66 0.79 1.00 7296 6471
##
## ~TrialNumber_c (Number of levels: 50)
## Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
## sd(Intercept) 0.09 0.03 0.04 0.15 1.00 1649 2746
##
## Regression Coefficients:
## Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
## Intercept 3.92 0.08 3.75 4.08 1.00 1736 3195
## CommAtt2M1 2.78 0.11 2.56 2.99 1.00 2927 3842
## CommAtt3M2 -3.95 0.13 -4.19 -3.70 1.00 2998 3937
## Modality1 -1.11 0.03 -1.17 -1.06 1.00 13750 6634
## Modality2 0.56 0.03 0.51 0.62 1.00 15623 6248
## Big5 1.15 0.03 1.10 1.21 1.00 10082 6655
## Familiarity 1.13 0.03 1.07 1.18 1.00 9422 6474
## Expressibility_z 0.39 0.01 0.36 0.41 1.00 16992 6093
##
## Further Distributional Parameters:
## Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
## sigma 0.52 0.01 0.50 0.55 1.00 11160 6372
## nu 2.96 0.16 2.67 3.28 1.00 13245 5813
##
## Draws were sampled using sampling(NUTS). For each parameter, Bulk_ESS
## and Tail_ESS are effective sample size measures, and Rhat is the potential
## scale reduction factor on split chains (at convergence, Rhat = 1).# the coefficients look ok again. Familiarity and Big5 even got closer to our simulated betasplot(h1.m5)
# all good
plot(conditional_effects(h1.m5), points = TRUE)
# the effects all go in good direction
pp_check(h1.m5, type = "dens_overlay")
# the fit seems somewhat better than with gaussian, but still negative values are there
pp_check(h1.m5, type = "error_scatter_avg")
# unchanged
h1.m5_R2## Estimate Est.Error Q2.5 Q97.5
## R2 0.8541097 0.002443239 0.84921 0.8587499# explained variance remains around 85%Model 6 - log-normal distribution
We have already seen - when plotting - that effort probably tends towards lognormal distribution. We will keep the model constant, just exchange the distribution. Priors are kept default
h1.m6 <- brm(Eff ~ 1 + CommAtt + Modality + Big5 + Familiarity + Expressibility_z + (1 + CommAtt || Participant) + (1 + CommAtt || Concept) + (1 || TrialNumber_c),
data = final_data,
family = lognormal(),
iter = 4000,
cores = 4)
# Add criterions for later diagnostics
h1.m6 <- add_criterion(h1.m6, criterion = c("loo", "waic"))
# Calculate also variance explained (R^2)
h1.m6_R2 <- bayes_R2(h1.m6)
# Save both as objects
saveRDS(h1.m6, here("09_Analysis_Modeling", "models", "h1.m6.rds"))
saveRDS(h1.m6_R2, here("09_Analysis_Modeling", "models", "h1.m6_R2.rds"))
beep(5)h1.m6 <- readRDS(here("09_Analysis_Modeling", "models", "h1.m6.rds"))
h1.m6_R2 <- readRDS(here("09_Analysis_Modeling", "models", "h1.m6_R2.rds"))
# Summary
summary(h1.m6)## Family: lognormal
## Links: mu = identity; sigma = identity
## Formula: Eff ~ 1 + CommAtt + Modality + Big5 + Familiarity + Expressibility_z + (1 + CommAtt || Participant) + (1 + CommAtt || Concept) + (1 || TrialNumber_c)
## Data: final_data (Number of observations: 5076)
## Draws: 4 chains, each with iter = 4000; warmup = 2000; thin = 1;
## total post-warmup draws = 8000
##
## Multilevel Hyperparameters:
## ~Concept (Number of levels: 21)
## Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
## sd(Intercept) 0.01 0.00 0.00 0.02 1.00 2298 2414
## sd(CommAtt2M1) 0.00 0.00 0.00 0.01 1.00 4803 3682
## sd(CommAtt3M2) 0.01 0.01 0.00 0.02 1.00 3687 3650
##
## ~Participant (Number of levels: 120)
## Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
## sd(Intercept) 0.05 0.00 0.04 0.06 1.00 2844 4508
## sd(CommAtt2M1) 0.00 0.00 0.00 0.01 1.00 4417 3655
## sd(CommAtt3M2) 0.01 0.01 0.00 0.03 1.00 2618 2888
##
## ~TrialNumber_c (Number of levels: 50)
## Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
## sd(Intercept) 0.00 0.00 0.00 0.01 1.00 3411 3578
##
## Regression Coefficients:
## Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
## Intercept 1.23 0.01 1.22 1.24 1.00 3120 4439
## CommAtt2M1 0.62 0.01 0.61 0.63 1.00 9634 6358
## CommAtt3M2 -1.10 0.01 -1.11 -1.08 1.00 9514 5588
## Modality1 -0.31 0.01 -0.33 -0.30 1.00 9097 6088
## Modality2 0.16 0.01 0.15 0.17 1.00 9223 6512
## Big5 0.32 0.01 0.31 0.34 1.00 2924 4558
## Familiarity 0.32 0.01 0.31 0.34 1.00 2297 3545
## Expressibility_z 0.12 0.00 0.12 0.13 1.00 10923 6710
##
## Further Distributional Parameters:
## Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
## sigma 0.17 0.00 0.17 0.18 1.00 11348 5998
##
## Draws were sampled using sampling(NUTS). For each parameter, Bulk_ESS
## and Tail_ESS are effective sample size measures, and Rhat is the potential
## scale reduction factor on split chains (at convergence, Rhat = 1).# Coefficients remain mostly unchanged Converting to original scale
Of course now we have all on lognormal scale so it is not so straightforwadly interpretable. This is how we can get to the values of an effort under certain value of a predictor
library(emmeans)
# Compute estimated marginal means on log-scale
em_h1.m6 <- emmeans(h1.m6, ~ CommAtt) #~ CommAtt
#Backtransform the post.beta values
em_h1.m6@post.beta <- exp(em_h1.m6@post.beta)
print(em_h1.m6)## CommAtt emmean lower.HPD upper.HPD
## 1 3.33 3.29 3.37
## 2 6.21 6.13 6.29
## 3 2.08 2.05 2.11
##
## Results are averaged over the levels of: Modality
## Point estimate displayed: median
## HPD interval probability: 0.95So we indeed see that effort in CommAtt2 increase but then decreases again for CommAtt3. We can also try to get our simulated betas
coeff1 <- 3.33*1.5
coeff2 <- 3.33*0.5
print(coeff1)## [1] 4.995print(coeff2) # which is very close to the mean values we see from the model## [1] 1.665We see that the effect estimated by the model is now stronger than in our simulated data - also maybe because the values are averaged over predictor modality
Now let’s try a different method to get the coefficients (code adapted from https://bruno.nicenboim.me/bayescogsci/ch-reg.html#sec-trial)
# Extract posterior samples
alpha_samples <- as_draws_df(h1.m6)$b_Intercept
beta_2_vs_1 <- as_draws_df(h1.m6)$b_CommAtt2M1
beta_3_vs_2 <- as_draws_df(h1.m6)$b_CommAtt3M2
# Compute expected values on the log scale
mu_1 <- alpha_samples # CommAtt 1
mu_2 <- alpha_samples + beta_2_vs_1 # CommAtt 2
mu_3 <- alpha_samples + beta_2_vs_1 + beta_3_vs_2 # CommAtt 3
# Transform to original scale
effect_1 <- exp(mu_1)
effect_2 <- exp(mu_2)
effect_3 <- exp(mu_3)
# Calculate contrasts on the original scale
effect_diff_2_vs_1 <- effect_2 - effect_1
effect_diff_3_vs_2 <- effect_3 - effect_2
effect_diff_3_vs_1 <- effect_3 - effect_1
# Summarize the effects
list(
mean_intercept = mean(effect_1),
mean_diff_2_vs_1 = c(mean = mean(effect_diff_2_vs_1), quantile(effect_diff_2_vs_1, c(0.025, 0.975))),
mean_diff_3_vs_2 = c(mean = mean(effect_diff_3_vs_2), quantile(effect_diff_3_vs_2, c(0.025, 0.975))),
mean_diff_3_vs_1 = c(mean = mean(effect_diff_3_vs_1), quantile(effect_diff_3_vs_1, c(0.025, 0.975)))
)## $mean_intercept
## [1] 3.419549
##
## $mean_diff_2_vs_1
## mean 2.5% 97.5%
## 2.950688 2.872812 3.029410
##
## $mean_diff_3_vs_2
## mean 2.5% 97.5%
## -4.240333 -4.321665 -4.160104
##
## $mean_diff_3_vs_1
## mean 2.5% 97.5%
## -1.289646 -1.320918 -1.257698Now we can use these coefficients to transform to our simulated betas
h1m6_coeff <- transform_attempt(3.419549, 2.950688, -4.240333)
print(h1m6_coeff)## b_attempt2 b_attempt3
## 1 1.862888 0.3343524Now this looks closer again to our betas.
This is for all predictors (except concept and participant)
# Extract posterior samples
posterior_samples <- as_draws_df(h1.m6)
alpha_samples <- posterior_samples$b_Intercept
# Create a list to store effects for each fixed factor
effect_list <- list()
# Helper function to calculate summary statistics
get_effect_summary <- function(effect_samples) {
mean_effect <- mean(effect_samples)
se_effect <- sd(effect_samples)
ci_effect <- quantile(effect_samples, c(0.025, 0.975))
post_prob <- mean(effect_samples > 0)
c(mean = mean_effect,
se = se_effect,
lower_ci = ci_effect[1],
upper_ci = ci_effect[2],
post_prob = post_prob)
}
# COMMATT (successive differences coding)
if ("b_CommAtt2M1" %in% colnames(posterior_samples) & "b_CommAtt3M2" %in% colnames(posterior_samples)) {
beta_2_vs_1 <- posterior_samples$b_CommAtt2M1
beta_3_vs_2 <- posterior_samples$b_CommAtt3M2
mu_1 <- alpha_samples
mu_2 <- alpha_samples + beta_2_vs_1
mu_3 <- alpha_samples + beta_2_vs_1 + beta_3_vs_2
effect_list$CommAtt <- rbind(
"commat 2 vs 1" = get_effect_summary(exp(mu_2) - exp(mu_1)),
"commat 3 vs 2" = get_effect_summary(exp(mu_3) - exp(mu_2)),
"commat 3 vs 1" = get_effect_summary(exp(mu_3) - exp(mu_1))
)
}
# MODALITY (sum contrasts scaled by 0.5)
if ("b_Modality1" %in% colnames(posterior_samples) & "b_Modality2" %in% colnames(posterior_samples)) {
beta_mod_1 <- posterior_samples$b_Modality1
beta_mod_2 <- posterior_samples$b_Modality2
mu_mod_1 <- alpha_samples + beta_mod_1
mu_mod_2 <- alpha_samples + beta_mod_2
mu_mod_3 <- alpha_samples - beta_mod_1 - beta_mod_2
effect_list$Modality <- rbind(
"mod 1 vs 2" = get_effect_summary(exp(mu_mod_1) - exp(mu_mod_2)),
"mod 1 vs 3" = get_effect_summary(exp(mu_mod_1) - exp(mu_mod_3)),
"mod 2 vs 3" = get_effect_summary(exp(mu_mod_2) - exp(mu_mod_3))
)
}
# BIG5 (continuous)
if ("b_Big5" %in% colnames(posterior_samples)) {
beta_big5 <- posterior_samples$b_Big5
effect_list$Big5 <- get_effect_summary(exp(alpha_samples + beta_big5) - exp(alpha_samples))
}
# FAMILIARITY (continuous)
if ("b_Familiarity" %in% colnames(posterior_samples)) {
beta_fam <- posterior_samples$b_Familiarity
effect_list$Familiarity <- get_effect_summary(exp(alpha_samples + beta_fam) - exp(alpha_samples))
}
# EXPRESSIBILITY_Z (continuous)
if ("b_Expressibility_z" %in% colnames(posterior_samples)) {
beta_expr <- posterior_samples$b_Expressibility_z
effect_list$Expressibility_z <- get_effect_summary(exp(alpha_samples + beta_expr) - exp(alpha_samples))
}
# TRIAL NUMBER (centered continuous)
if ("b_TrialNumber_c" %in% colnames(posterior_samples)) {
beta_trial <- posterior_samples$b_TrialNumber_c
effect_list$TrialNumber_c <- get_effect_summary(exp(alpha_samples + beta_trial) - exp(alpha_samples))
}
# Convert to a nicely formatted data frame
effect_df <- do.call(rbind, effect_list)
# View effects
effect_df## mean se lower_ci.2.5% upper_ci.97.5% post_prob
## commat 2 vs 1 2.95068779 0.04052085 2.87281233 3.0294096 1.00000
## commat 3 vs 2 -4.24033346 0.04140812 -4.32166524 -4.1601045 0.00000
## commat 3 vs 1 -1.28964568 0.01620638 -1.32091801 -1.2576981 0.00000
## mod 1 vs 2 -1.52303929 0.04030161 -1.60299691 -1.4435554 0.00000
## mod 1 vs 3 -1.49131528 0.04098655 -1.57288685 -1.4107477 0.00000
## mod 2 vs 3 0.03172401 0.04858764 -0.06380552 0.1252410 0.74175
## Big5 1.31309529 0.04264827 1.23007768 1.3970338 1.00000
## Familiarity 1.30972571 0.04408809 1.22420946 1.3983930 1.00000
## Expressibility_z 0.44448503 0.01197175 0.42135330 0.4684496 1.00000So here we also see the negative effect of combined modality (mod1)
plot(h1.m6)
# all good
plot(conditional_effects(h1.m6), points = TRUE)
# the effects all go in good direction
pp_check(h1.m6, type = "dens_overlay")
# we got rid of negative predictions, and it looks very good
pp_check(h1.m6, type = "error_scatter_avg")
# very nice blob
h1.m6_R2## Estimate Est.Error Q2.5 Q97.5
## R2 0.8889641 0.00137284 0.8861744 0.8915303# explained variance increases to 88%Model 6.1 - with correlation
Since we now significantly improved the model performance, let’s try once again the correlation between slope and intercept
h1.m6c <- brm(Eff ~ 1 + CommAtt + Modality + Big5 + Familiarity + Expressibility_z + (1 + CommAtt | Participant) + (1 + CommAtt | Concept) + (1 | TrialNumber_c),
data = final_data,
family = lognormal(),
iter = 4000,
cores = 4)
# Add criterions for later diagnostics
h1.m6c <- add_criterion(h1.m6c, criterion = c("loo", "waic"))
# Calculate also variance explained (R^2)
h1.m6c_R2 <- bayes_R2(h1.m6c)
# Save both as objects
saveRDS(h1.m6c, here("09_Analysis_Modeling", "models", "h1.m6c.rds"))
saveRDS(h1.m6c_R2, here("09_Analysis_Modeling", "models", "h1.m6c_R2.rds"))
beep(5)h1.m6c <- readRDS(here("09_Analysis_Modeling", "models", "h1.m6c.rds"))
h1.m6c_R2 <- readRDS(here("09_Analysis_Modeling", "models", "h1.m6c_R2.rds"))
# Summary
summary(h1.m6c)## Family: lognormal
## Links: mu = identity; sigma = identity
## Formula: Eff ~ 1 + CommAtt + Modality + Big5 + Familiarity + Expressibility_z + (1 + CommAtt | Participant) + (1 + CommAtt | Concept) + (1 | TrialNumber_c)
## Data: final_data (Number of observations: 5076)
## Draws: 4 chains, each with iter = 4000; warmup = 2000; thin = 1;
## total post-warmup draws = 8000
##
## Multilevel Hyperparameters:
## ~Concept (Number of levels: 21)
## Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS
## sd(Intercept) 0.01 0.00 0.00 0.02 1.00 1744
## sd(CommAtt2M1) 0.01 0.00 0.00 0.02 1.00 5613
## sd(CommAtt3M2) 0.01 0.01 0.00 0.03 1.00 3784
## cor(Intercept,CommAtt2M1) 0.28 0.49 -0.78 0.95 1.00 11362
## cor(Intercept,CommAtt3M2) 0.33 0.47 -0.75 0.95 1.00 7628
## cor(CommAtt2M1,CommAtt3M2) -0.05 0.50 -0.89 0.87 1.00 8530
## Tail_ESS
## sd(Intercept) 1448
## sd(CommAtt2M1) 5243
## sd(CommAtt3M2) 4533
## cor(Intercept,CommAtt2M1) 5945
## cor(Intercept,CommAtt3M2) 5887
## cor(CommAtt2M1,CommAtt3M2) 6522
##
## ~Participant (Number of levels: 120)
## Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS
## sd(Intercept) 0.05 0.00 0.04 0.06 1.00 3640
## sd(CommAtt2M1) 0.01 0.00 0.00 0.02 1.00 5290
## sd(CommAtt3M2) 0.01 0.01 0.00 0.03 1.00 3474
## cor(Intercept,CommAtt2M1) 0.35 0.46 -0.73 0.95 1.00 10725
## cor(Intercept,CommAtt3M2) -0.40 0.40 -0.94 0.59 1.00 9703
## cor(CommAtt2M1,CommAtt3M2) -0.27 0.49 -0.94 0.78 1.00 5410
## Tail_ESS
## sd(Intercept) 5384
## sd(CommAtt2M1) 4470
## sd(CommAtt3M2) 4069
## cor(Intercept,CommAtt2M1) 5934
## cor(Intercept,CommAtt3M2) 5524
## cor(CommAtt2M1,CommAtt3M2) 6511
##
## ~TrialNumber_c (Number of levels: 50)
## Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
## sd(Intercept) 0.00 0.00 0.00 0.01 1.00 4362 4575
##
## Regression Coefficients:
## Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
## Intercept 1.23 0.01 1.22 1.24 1.00 5205 5979
## CommAtt2M1 0.62 0.01 0.61 0.63 1.00 14611 6124
## CommAtt3M2 -1.10 0.01 -1.11 -1.08 1.00 12219 6336
## Modality1 -0.32 0.01 -0.33 -0.30 1.00 14134 6778
## Modality2 0.16 0.01 0.15 0.18 1.00 15159 6553
## Big5 0.32 0.01 0.31 0.34 1.00 4826 6410
## Familiarity 0.32 0.01 0.30 0.34 1.00 4174 5246
## Expressibility_z 0.12 0.00 0.12 0.13 1.00 13884 6876
##
## Further Distributional Parameters:
## Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
## sigma 0.17 0.00 0.17 0.18 1.00 17858 5344
##
## Draws were sampled using sampling(NUTS). For each parameter, Bulk_ESS
## and Tail_ESS are effective sample size measures, and Rhat is the potential
## scale reduction factor on split chains (at convergence, Rhat = 1).# Coefficients remain mostly unchanged plot(h1.m6c)
# now correlation does not seem to generate problems
plot(conditional_effects(h1.m6c), points = TRUE)
# the effects all go in good direction
pp_check(h1.m6c, type = "dens_overlay")
# ppcheck good
pp_check(h1.m6c, type = "error_scatter_avg")
# nice blob
h1.m6c_R2## Estimate Est.Error Q2.5 Q97.5
## R2 0.8890392 0.001393478 0.8861745 0.8916361# explained variance remains around 88%Diagnostics I
Let’s now check which model has the best diagnostics
Rhat
# Extract R-hat values for each model
rhat_list <- lapply(model_list, function(model) {
rhat_values <- rhat(model)
data.frame(model = deparse(substitute(model)),
max_rhat = max(rhat_values),
min_rhat = min(rhat_values))
})
# Combine and inspect
do.call(rbind, rhat_list)## model max_rhat min_rhat
## 1 X[[i]] 1.002349 0.9996758
## 2 X[[i]] 1.011516 0.9997215
## 3 X[[i]] 1.005103 0.9997167
## 4 X[[i]] 1.003503 0.9996573
## 5 X[[i]] 1.004119 0.9996640
## 6 X[[i]] 1.003712 0.9996386
## 7 X[[i]] 1.002459 0.9995773
## 8 X[[i]] 1.002605 0.9996426All models seems actually ok in terms of Rhat values except model 2 (h1.m2)
ESS
Effective sample size tells how many independent samples the model has effectively drawn from the PD. Low ESS suggests autocorrelation (i.e., sample explores one part of posterior), while high ESS means good mix
# Extract n_eff values for each model
neff_ratio_list <- lapply(model_list, function(model) {
neff_values <- neff_ratio(model) # Here we calculate ratio (not the raw number of effective samples)
data.frame(model = deparse(substitute(model)),
min_neff = min(neff_values),
max_neff = max(neff_values),
mean_neff = mean(neff_values))
})
# Combine and inspect
do.call(rbind, neff_ratio_list)## model min_neff max_neff mean_neff
## 1 X[[i]] 0.23283181 0.8832823 0.6862927
## 2 X[[i]] 0.04396011 0.9020581 0.6204512
## 3 X[[i]] 0.14481431 0.8834718 0.7504325
## 4 X[[i]] 0.22967606 0.9020381 0.7709286
## 5 X[[i]] 0.11769241 0.9065446 0.7505601
## 6 X[[i]] 0.20616401 0.9322803 0.7457049
## 7 X[[i]] 0.28072060 0.8984722 0.7626900
## 8 X[[i]] 0.18101711 0.9370772 0.7924409So the highest ratio have model h1.m6c (lognormal with correlation) but in fact they are all quite comparable. Let’s loot at 3 highest
effective_sample(h1.m6c) # this one indeed seems the best as it has all ESS around 10k## Parameter ESS
## 1 b_Intercept 5193
## 2 b_CommAtt2M1 14438
## 3 b_CommAtt3M2 12340
## 4 b_Modality1 14518
## 5 b_Modality2 14927
## 6 b_Big5 4824
## 7 b_Familiarity 4162
## 8 b_Expressibility_z 14233effective_sample(h1.m6)## Parameter ESS
## 1 b_Intercept 3104
## 2 b_CommAtt2M1 9618
## 3 b_CommAtt3M2 9509
## 4 b_Modality1 9057
## 5 b_Modality2 9192
## 6 b_Big5 2906
## 7 b_Familiarity 2280
## 8 b_Expressibility_z 10811effective_sample(h1.m3p) ## Parameter ESS
## 1 b_Intercept 8533
## 2 b_CommAtt2M1 2795
## 3 b_CommAtt3M2 2385
## 4 b_Modality1 13744
## 5 b_Modality2 13144
## 6 b_Big5 12483
## 7 b_Familiarity 13139
## 8 b_Expressibility_z 15614So there are some considerable differences for different coefficients
LOO & WAIC
Leave-one-out (loo) validation
l <- loo_compare(h1.m1, h1.m2, h1.m3, h1.m3p, h1.m4, h1.m5, h1.m6, h1.m6c, criterion = "loo")
print(l, simplify = F)## elpd_diff se_diff elpd_loo se_elpd_loo p_loo se_p_loo looic se_looic
## h1.m6c 0.0 0.0 -5070.5 69.3 126.5 3.5 10141.0 138.5
## h1.m6 -0.2 1.5 -5070.7 69.2 120.4 3.3 10141.4 138.5
## h1.m5 -904.2 59.2 -5974.6 74.1 418.5 5.3 11949.3 148.1
## h1.m2 -1143.5 67.5 -6214.0 78.4 285.6 10.7 12428.0 156.8
## h1.m3 -1221.1 66.5 -6291.6 75.5 261.3 7.7 12583.2 150.9
## h1.m3p -1222.5 66.9 -6293.0 76.3 281.6 8.3 12586.0 152.5
## h1.m4 -1224.7 67.4 -6295.2 77.0 281.3 8.4 12590.3 154.0
## h1.m1 -1681.0 72.2 -6751.5 74.9 25.3 0.8 13503.0 149.9elpd_loo: This is the expected log pointwise predictive density for LOO. Higher values indicate a better fit to the data.
se_elpd_loo: The standard error of the elpd_loo, representing uncertainty in the model’s predictive fit according to LOO.
looic: The LOO Information Criterion, which is similar to waic but based on leave-one-out cross-validation. Lower values are better.
p_loo: The effective number of parameters according to LOO, indicating the model’s complexity.
se_p_loo: The standard error of p_loo, representing uncertainty around the effective number of parameters.
So lognormal seems the best.
Information criterion (WAIC)
w <- loo_compare(h1.m1, h1.m2, h1.m3, h1.m3p, h1.m4, h1.m5, h1.m6, h1.m6c, criterion = "waic")
print(w, simplify = F)## elpd_diff se_diff elpd_waic se_elpd_waic p_waic se_p_waic waic
## h1.m6c 0.0 0.0 -5070.3 69.3 126.3 3.5 10140.5
## h1.m6 -0.3 1.5 -5070.5 69.2 120.2 3.3 10141.0
## h1.m5 -903.7 59.2 -5974.0 74.1 417.8 5.3 11948.0
## h1.m2 -1142.2 67.5 -6212.5 78.3 284.1 10.6 12424.9
## h1.m3 -1220.4 66.5 -6290.7 75.4 260.3 7.6 12581.3
## h1.m3p -1221.6 66.9 -6291.8 76.2 280.4 8.2 12583.7
## h1.m4 -1223.8 67.3 -6294.1 77.0 280.3 8.4 12588.2
## h1.m1 -1681.2 72.2 -6751.5 74.9 25.3 0.8 13502.9
## se_waic
## h1.m6c 138.5
## h1.m6 138.5
## h1.m5 148.1
## h1.m2 156.6
## h1.m3 150.9
## h1.m3p 152.5
## h1.m4 153.9
## h1.m1 149.9# see Solomon Kurz
cbind(waic_diff = w[,1] * -2,
se = w[,2] * 2)## waic_diff se
## h1.m6c 0.0000000 0.000000
## h1.m6 0.5289111 2.917133
## h1.m5 1807.4580148 118.471557
## h1.m2 2284.3843630 134.987550
## h1.m3 2440.8062251 132.903242
## h1.m3p 2443.1522534 133.864984
## h1.m4 2447.6767321 134.678147
## h1.m1 3362.3866702 144.475730elpd_waic (expected log pointwise predictive density for WAIC): This represents the model’s predictive fit to the data. Higher values indicate a better fit.
se_elpd_waic (standard error of elpd_waic): Measures uncertainty around the elpd_waic estimate.
waic: The Widely Applicable Information Criterion, a measure of model fit where lower values indicate a better fit.
se_waic (standard error of WAIC): Uncertainty around the WAIC estimate.
elpd_diff: The difference in the elpd_waic between the model in question and the baseline model (fit_eff_2, which has elpd_diff of 0). A negative value indicates that the model fits worse than fit_eff_2.
se_diff: The standard error of the elpd_diff, indicating how much uncertainty there is in the difference in predictive performance.
p_waic: The number of effective parameters in the model (related to model complexity). Lower values indicate simpler models, and higher values suggest more complexity.
Plot the comparison
model_weights(h1.m1, h1.m2, h1.m3, h1.m3p, h1.m4, h1.m5, h1.m6, h1.m6c, weights = "waic") %>%
round(digits = 2)## h1.m1 h1.m2 h1.m3 h1.m3p h1.m4 h1.m5 h1.m6 h1.m6c
## 0.00 0.00 0.00 0.00 0.00 0.00 0.43 0.57So as ppcheck already suggested, lognormal model indeed seem to have the most predictive power. For this particular (synthetic) data, we will now proceed with model h1.m6c
We will first add some mildly informative priors, and then we also try to add some interaction terms and do a comparison once again.
Model 7 - lognormal with priors
Let’s first check what priors have been selected as default for h1.m6c
# Print priors
prior_summary(h1.m6c)## prior class coef group resp dpar
## (flat) b
## (flat) b Big5
## (flat) b CommAtt2M1
## (flat) b CommAtt3M2
## (flat) b Expressibility_z
## (flat) b Familiarity
## (flat) b Modality1
## (flat) b Modality2
## student_t(3, 1.3, 2.5) Intercept
## lkj_corr_cholesky(1) L
## lkj_corr_cholesky(1) L Concept
## lkj_corr_cholesky(1) L Participant
## student_t(3, 0, 2.5) sd
## student_t(3, 0, 2.5) sd Concept
## student_t(3, 0, 2.5) sd CommAtt2M1 Concept
## student_t(3, 0, 2.5) sd CommAtt3M2 Concept
## student_t(3, 0, 2.5) sd Intercept Concept
## student_t(3, 0, 2.5) sd Participant
## student_t(3, 0, 2.5) sd CommAtt2M1 Participant
## student_t(3, 0, 2.5) sd CommAtt3M2 Participant
## student_t(3, 0, 2.5) sd Intercept Participant
## student_t(3, 0, 2.5) sd TrialNumber_c
## student_t(3, 0, 2.5) sd Intercept TrialNumber_c
## student_t(3, 0, 2.5) sigma
## nlpar lb ub source
## default
## (vectorized)
## (vectorized)
## (vectorized)
## (vectorized)
## (vectorized)
## (vectorized)
## (vectorized)
## default
## default
## (vectorized)
## (vectorized)
## 0 default
## 0 (vectorized)
## 0 (vectorized)
## 0 (vectorized)
## 0 (vectorized)
## 0 (vectorized)
## 0 (vectorized)
## 0 (vectorized)
## 0 (vectorized)
## 0 (vectorized)
## 0 (vectorized)
## 0 defaultOk, we can keep all defaulted ones, but we do not need to leave flat priors for the beta coefficients as we do have some assumptions/expectations
Let’s reuse priors we have already used for h1.m3p
priors_h1m7 <- c(
set_prior("normal(2.5, 0.5)", class = "Intercept", lb=0),
set_prior("normal(0,0.50)", class = "b", coef = "CommAtt2M1"),
set_prior("normal(0,0.50)", class = "b", coef = "CommAtt3M2"),
set_prior("normal(0,0.25)", class = "b", coef = "Modality1"),
set_prior("normal(0,0.25)", class = "b", coef = "Modality2"),
set_prior("normal(0,0.25)", class = "b", coef = "Big5"),
set_prior("normal(0,0.25)", class = "b", coef = "Familiarity"),
set_prior("normal(0,0.25)", class = "b", coef = "Expressibility_z")
)
# The rest we will leave default (and check afterwards)
h1.m7 <- brm(Eff ~ 1 + CommAtt + Modality + Big5 + Familiarity + Expressibility_z + (1 + CommAtt | Participant) + (1 + CommAtt | Concept) + (1 | TrialNumber_c),
data = final_data,
family = lognormal(),
prior = priors_h1m7,
iter = 4000,
cores = 4)
# Add criterions for later diagnostics
h1.m7 <- add_criterion(h1.m7, criterion = c("loo", "waic"))
# Calculate also variance explained (R^2)
h1.m7_R2 <- bayes_R2(h1.m7)
# Save both as objects
saveRDS(h1.m7, here("09_Analysis_Modeling", "models", "h1.m7.rds"))
saveRDS(h1.m7_R2, here("09_Analysis_Modeling", "models", "h1.m7_R2.rds"))
beep(5)h1.m7 <- readRDS(here("09_Analysis_Modeling", "models", "h1.m7.rds"))
h1.m7_R2 <- readRDS(here("09_Analysis_Modeling", "models", "h1.m7_R2.rds"))
# Summary
summary(h1.m7)## Family: lognormal
## Links: mu = identity; sigma = identity
## Formula: Eff ~ 1 + CommAtt + Modality + Big5 + Familiarity + Expressibility_z + (1 + CommAtt | Participant) + (1 + CommAtt | Concept) + (1 | TrialNumber_c)
## Data: final_data (Number of observations: 5076)
## Draws: 4 chains, each with iter = 4000; warmup = 2000; thin = 1;
## total post-warmup draws = 8000
##
## Multilevel Hyperparameters:
## ~Concept (Number of levels: 21)
## Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS
## sd(Intercept) 0.01 0.00 0.00 0.02 1.00 2406
## sd(CommAtt2M1) 0.01 0.00 0.00 0.02 1.00 5336
## sd(CommAtt3M2) 0.01 0.01 0.00 0.03 1.00 3663
## cor(Intercept,CommAtt2M1) 0.27 0.49 -0.79 0.95 1.00 9049
## cor(Intercept,CommAtt3M2) 0.34 0.47 -0.74 0.95 1.00 7153
## cor(CommAtt2M1,CommAtt3M2) -0.05 0.50 -0.89 0.88 1.00 8030
## Tail_ESS
## sd(Intercept) 2582
## sd(CommAtt2M1) 4860
## sd(CommAtt3M2) 4212
## cor(Intercept,CommAtt2M1) 5794
## cor(Intercept,CommAtt3M2) 6108
## cor(CommAtt2M1,CommAtt3M2) 5937
##
## ~Participant (Number of levels: 120)
## Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS
## sd(Intercept) 0.05 0.00 0.04 0.06 1.00 3562
## sd(CommAtt2M1) 0.01 0.00 0.00 0.02 1.00 5277
## sd(CommAtt3M2) 0.01 0.01 0.00 0.03 1.00 3118
## cor(Intercept,CommAtt2M1) 0.34 0.46 -0.72 0.95 1.00 11493
## cor(Intercept,CommAtt3M2) -0.40 0.40 -0.94 0.57 1.00 7854
## cor(CommAtt2M1,CommAtt3M2) -0.26 0.48 -0.94 0.77 1.00 5771
## Tail_ESS
## sd(Intercept) 5189
## sd(CommAtt2M1) 4743
## sd(CommAtt3M2) 3237
## cor(Intercept,CommAtt2M1) 6204
## cor(Intercept,CommAtt3M2) 5447
## cor(CommAtt2M1,CommAtt3M2) 6837
##
## ~TrialNumber_c (Number of levels: 50)
## Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
## sd(Intercept) 0.00 0.00 0.00 0.01 1.00 4640 4109
##
## Regression Coefficients:
## Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
## Intercept 1.23 0.01 1.22 1.24 1.00 5166 5351
## CommAtt2M1 0.62 0.01 0.61 0.63 1.00 12440 6165
## CommAtt3M2 -1.10 0.01 -1.11 -1.08 1.00 10543 6212
## Modality1 -0.31 0.01 -0.33 -0.30 1.00 13106 5622
## Modality2 0.16 0.01 0.15 0.17 1.00 13551 6493
## Big5 0.32 0.01 0.31 0.34 1.00 4012 4552
## Familiarity 0.32 0.01 0.30 0.34 1.00 4337 5552
## Expressibility_z 0.12 0.00 0.12 0.13 1.00 11948 6505
##
## Further Distributional Parameters:
## Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
## sigma 0.17 0.00 0.17 0.18 1.00 15087 5361
##
## Draws were sampled using sampling(NUTS). For each parameter, Bulk_ESS
## and Tail_ESS are effective sample size measures, and Rhat is the potential
## scale reduction factor on split chains (at convergence, Rhat = 1).# Coefficients remain mostly unchanged plot(h1.m7)
# all good
plot(conditional_effects(h1.m7), points = TRUE)
# the effects all go in good direction
pp_check(h1.m7, type = "dens_overlay")
# nice
pp_check(h1.m7, type = "error_scatter_avg")
# unchanged
h1.m7_R2## Estimate Est.Error Q2.5 Q97.5
## R2 0.8889672 0.001399296 0.8861732 0.8916425# explained variance remains around 88%Let’s now also check whether the priors make sense
Model 8 - adding interactions
Possible interactions include:
- CommAtt x Modality - Effort could increase differently across modalities, depending on whether concept is guessed on the first, second or third attempt. E.g., gesture might require more effort on initial attempt, but vocal require more effort in repeated attempt
Note: Interesting, it would help disentagle the benefit of multimodality over other - maybe more effort overal but less attempts. However, we already model effect of modality on effort, so maybe this is not top priority
- CommAtt x Expressibility - Higher expressibility should moderate the effect of repeated attempts, such that the increase in effort with each additional attempt is smaller (or bigger?) for more expressible concepts.
Note: Not priority I would say
- Modality × Expressibility_z - The influence of expressibility on effort could be modality-specific — perhaps effort increases less with expressibility in the combined modality.
Note: Not priority (especially since expressibiliy has already modality encoded)
- Familiarity x CommAtt - More familiar partners may guess faster (fewer attempts) and require less effort, but this effect could diminish over multiple attempts.
Note: Does not seem to be priority
- Big5 × Modality or Big5 × CommAtt - More open/extraverted participants might maintain higher effort over attempts, or adjust more dynamically depending on the communicative channel.
Note: Not priority, but it is interesting if we want to tap more into the interindividual variability
Let’s try CommAtt x Modality and Big5 x CommAtt
h1.m8 <- brm(Eff ~ 1 + CommAtt * Modality + Big5 * CommAtt + Familiarity + Expressibility_z + (1 + CommAtt | Participant) + (1 + CommAtt | Concept) + (1 | TrialNumber_c),
data = final_data,
family = lognormal(),
prior = priors_h1m7, # we keep the priors from previous model
iter = 4000,
cores = 4)
# Add criterions for later diagnostics
h1.m8 <- add_criterion(h1.m8, criterion = c("loo", "waic"))
# Calculate also variance explained (R^2)
h1.m8_R2 <- bayes_R2(h1.m8)
# Save both as objects
saveRDS(h1.m8, here("09_Analysis_Modeling", "models", "h1.m8.rds"))
saveRDS(h1.m8_R2, here("09_Analysis_Modeling", "models", "h1.m8_R2.rds"))
beep(5)h1.m8 <- readRDS(here("09_Analysis_Modeling", "models", "h1.m8.rds"))
h1.m8_R2 <- readRDS(here("09_Analysis_Modeling", "models", "h1.m8_R2.rds"))
# Summary
summary(h1.m8)## Family: lognormal
## Links: mu = identity; sigma = identity
## Formula: Eff ~ 1 + CommAtt * Modality + Big5 * CommAtt + Familiarity + Expressibility_z + (1 + CommAtt | Participant) + (1 + CommAtt | Concept) + (1 | TrialNumber_c)
## Data: final_data (Number of observations: 5076)
## Draws: 4 chains, each with iter = 4000; warmup = 2000; thin = 1;
## total post-warmup draws = 8000
##
## Multilevel Hyperparameters:
## ~Concept (Number of levels: 21)
## Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS
## sd(Intercept) 0.01 0.00 0.00 0.02 1.00 2549
## sd(CommAtt2M1) 0.01 0.00 0.00 0.02 1.00 5836
## sd(CommAtt3M2) 0.01 0.01 0.00 0.03 1.00 4322
## cor(Intercept,CommAtt2M1) 0.29 0.48 -0.76 0.94 1.00 10754
## cor(Intercept,CommAtt3M2) 0.36 0.46 -0.70 0.95 1.00 8535
## cor(CommAtt2M1,CommAtt3M2) -0.04 0.50 -0.89 0.86 1.00 7971
## Tail_ESS
## sd(Intercept) 2905
## sd(CommAtt2M1) 4937
## sd(CommAtt3M2) 4495
## cor(Intercept,CommAtt2M1) 5912
## cor(Intercept,CommAtt3M2) 6615
## cor(CommAtt2M1,CommAtt3M2) 6944
##
## ~Participant (Number of levels: 120)
## Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS
## sd(Intercept) 0.05 0.00 0.04 0.06 1.00 3682
## sd(CommAtt2M1) 0.01 0.00 0.00 0.02 1.00 5068
## sd(CommAtt3M2) 0.01 0.01 0.00 0.03 1.00 3392
## cor(Intercept,CommAtt2M1) 0.34 0.45 -0.71 0.94 1.00 12364
## cor(Intercept,CommAtt3M2) -0.41 0.39 -0.94 0.57 1.00 9996
## cor(CommAtt2M1,CommAtt3M2) -0.26 0.49 -0.95 0.78 1.00 5862
## Tail_ESS
## sd(Intercept) 5433
## sd(CommAtt2M1) 5479
## sd(CommAtt3M2) 3982
## cor(Intercept,CommAtt2M1) 5945
## cor(Intercept,CommAtt3M2) 6038
## cor(CommAtt2M1,CommAtt3M2) 6673
##
## ~TrialNumber_c (Number of levels: 50)
## Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
## sd(Intercept) 0.00 0.00 0.00 0.01 1.00 4900 4106
##
## Regression Coefficients:
## Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS
## Intercept 1.23 0.01 1.22 1.24 1.00 5217
## CommAtt2M1 0.62 0.01 0.61 0.63 1.00 14494
## CommAtt3M2 -1.10 0.01 -1.11 -1.08 1.00 13959
## Modality1 -0.31 0.01 -0.33 -0.30 1.00 12978
## Modality2 0.16 0.01 0.15 0.18 1.00 13925
## Big5 0.32 0.01 0.30 0.34 1.00 4506
## Familiarity 0.32 0.01 0.30 0.34 1.00 3885
## Expressibility_z 0.12 0.00 0.12 0.13 1.00 12340
## CommAtt2M1:Modality1 0.01 0.02 -0.02 0.04 1.00 13065
## CommAtt3M2:Modality1 0.01 0.02 -0.03 0.05 1.00 12851
## CommAtt2M1:Modality2 -0.01 0.02 -0.04 0.02 1.00 13667
## CommAtt3M2:Modality2 0.02 0.02 -0.02 0.06 1.00 12109
## CommAtt2M1:Big5 -0.00 0.01 -0.02 0.02 1.00 14725
## CommAtt3M2:Big5 -0.01 0.01 -0.04 0.01 1.00 14109
## Tail_ESS
## Intercept 6099
## CommAtt2M1 7032
## CommAtt3M2 6248
## Modality1 6662
## Modality2 6627
## Big5 5386
## Familiarity 5579
## Expressibility_z 6403
## CommAtt2M1:Modality1 6974
## CommAtt3M2:Modality1 6884
## CommAtt2M1:Modality2 6960
## CommAtt3M2:Modality2 6404
## CommAtt2M1:Big5 6229
## CommAtt3M2:Big5 6504
##
## Further Distributional Parameters:
## Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
## sigma 0.17 0.00 0.17 0.18 1.00 17574 5831
##
## Draws were sampled using sampling(NUTS). For each parameter, Bulk_ESS
## and Tail_ESS are effective sample size measures, and Rhat is the potential
## scale reduction factor on split chains (at convergence, Rhat = 1).# Coefficients remain mostly unchanged
# Since we did not really focused on the interactions during the simulation, we also don't have strong expectations here. But for the real data, there is a good reason to expect some meaningful valuesplot(h1.m8)
# all good
plot(conditional_effects(h1.m8), points = TRUE)
# the effects all go in good direction
# we can see here that combined modality remains moderated across all commatts
# and also big5 seems to matter the most in the second attempt
pp_check(h1.m8, type = "dens_overlay")
# the problem with predicting negative values remains
pp_check(h1.m8, type = "error_scatter_avg")
# unchanged
h1.m8_R2## Estimate Est.Error Q2.5 Q97.5
## R2 0.8889759 0.001459646 0.8859485 0.8916926# explained variance remains around 87%Diagnostics II
ESS
effective_sample(h1.m6) ## Parameter ESS
## 1 b_Intercept 3104
## 2 b_CommAtt2M1 9618
## 3 b_CommAtt3M2 9509
## 4 b_Modality1 9057
## 5 b_Modality2 9192
## 6 b_Big5 2906
## 7 b_Familiarity 2280
## 8 b_Expressibility_z 10811effective_sample(h1.m6c)## Parameter ESS
## 1 b_Intercept 5193
## 2 b_CommAtt2M1 14438
## 3 b_CommAtt3M2 12340
## 4 b_Modality1 14518
## 5 b_Modality2 14927
## 6 b_Big5 4824
## 7 b_Familiarity 4162
## 8 b_Expressibility_z 14233effective_sample(h1.m7) ## Parameter ESS
## 1 b_Intercept 5149
## 2 b_CommAtt2M1 12363
## 3 b_CommAtt3M2 10483
## 4 b_Modality1 13095
## 5 b_Modality2 13808
## 6 b_Big5 3976
## 7 b_Familiarity 4298
## 8 b_Expressibility_z 11917effective_sample(h1.m8) ## Parameter ESS
## 1 b_Intercept 5221
## 2 b_CommAtt2M1 14165
## 3 b_CommAtt3M2 13961
## 4 b_Modality1 12872
## 5 b_Modality2 13799
## 6 b_Big5 4484
## 7 b_Familiarity 3878
## 8 b_Expressibility_z 12296
## 9 b_CommAtt2M1:Modality1 13103
## 10 b_CommAtt3M2:Modality1 12880
## 11 b_CommAtt2M1:Modality2 13651
## 12 b_CommAtt3M2:Modality2 12303
## 13 b_CommAtt2M1:Big5 15040
## 14 b_CommAtt3M2:Big5 13963# Now h1m6 looks as the weakest, while h1m8 looks much better - but ESS for Intercept is for some reason still quite lowLOO & WAIC
Add the criteria to each model
l <- loo_compare(h1.m6, h1.m6c, h1.m7, h1.m8, criterion = "loo")
print(l, simplify = F)## elpd_diff se_diff elpd_loo se_elpd_loo p_loo se_p_loo looic se_looic
## h1.m6c 0.0 0.0 -5070.5 69.3 126.5 3.5 10141.0 138.5
## h1.m6 -0.2 1.5 -5070.7 69.2 120.4 3.3 10141.4 138.5
## h1.m7 -0.3 0.2 -5070.8 69.3 126.8 3.5 10141.5 138.6
## h1.m8 -4.4 1.9 -5074.8 69.2 132.9 3.7 10149.7 138.5w <- loo_compare(h1.m6, h1.m6c, h1.m7, h1.m8, criterion = "waic")
print(w, simplify = F)## elpd_diff se_diff elpd_waic se_elpd_waic p_waic se_p_waic waic
## h1.m6c 0.0 0.0 -5070.3 69.3 126.3 3.5 10140.5
## h1.m6 -0.3 1.5 -5070.5 69.2 120.2 3.3 10141.0
## h1.m7 -0.3 0.2 -5070.5 69.3 126.5 3.5 10141.1
## h1.m8 -4.4 1.9 -5074.6 69.2 132.7 3.7 10149.2
## se_waic
## h1.m6c 138.5
## h1.m6 138.5
## h1.m7 138.6
## h1.m8 138.5# see Solomon Kurz
cbind(waic_diff = w[,1] * -2,
se = w[,2] * 2)## waic_diff se
## h1.m6c 0.0000000 0.0000000
## h1.m6 0.5289111 2.9171328
## h1.m7 0.5582963 0.3803033
## h1.m8 8.7265395 3.7869214So in terms of WAIC&LOO, interactions do not really add predictive power. This might be specific for the synthetic data, as we did not explicitly focused on the interaction coefficients there. At the same time, the difference of h1.m8 from h1.m6c is not so significant. For now, we stop here and with the real data, we will proceed with the same evaluation to pick up the best model
H2: Data wrangling
Now we can need to account also for our second hypothesis, namely
H2: A higher degree of misunderstanding will require a performer to engage in more effortful correction.
The difficulty with our data structure is that PrevAn is a variable that has values only for second and third correction (i.e., there is no previous answer for the first performance). If we left it unchange, the model will eventually get rid of all NA values, which means we will loose the relationship between effort in second and first communicative attempt. To avoid this, we will create a new variable which we call Effort Change Ratio. We will simply calculate the ratio of change from effort in communicative attempt x to communicative attempt x+1. Like that, we will still get rid of communicative attempt 1, but the ratio that will be associated with CommAtt==2 already encapsulates the relationship towards this attempt.
# Calculate Effort Change Ratio
final_data_2 <- final_data %>%
group_by(Participant, Concept) %>%
mutate(
Effort_1 = Eff[CommAtt == 1][1], # Effort for attempt 1
Effort_2 = Eff[CommAtt == 2][1], # Effort for attempt 2
Effort_3 = Eff[CommAtt == 3][1], # Effort for attempt 3
Effort_Change_Ratio_1_to_2 = case_when(
CommAtt == 2 & !is.na(Effort_1) ~ Eff / Effort_1,
TRUE ~ NA_real_
),
Effort_Change_Ratio_2_to_3 = case_when(
CommAtt == 3 & !is.na(Effort_2) ~ Eff / Effort_2,
TRUE ~ NA_real_
)
) %>%
ungroup()
# Collide ratios into single col
final_data_2 <- final_data_2 %>%
mutate(
Effort_Change_Ratio = coalesce(Effort_Change_Ratio_1_to_2, Effort_Change_Ratio_2_to_3)
)
# Remove unnecessary columns
final_data_2 <- subset(final_data_2, select = -c(Effort_1, Effort_2, Effort_3, Effort_Change_Ratio_1_to_2, Effort_Change_Ratio_2_to_3))
# View the result
head(final_data_2)## # A tibble: 6 × 11
## Participant Concept Modality Big5 Familiarity Expressibility CommAtt Eff
## <dbl> <dbl> <chr> <dbl> <dbl> <dbl> <dbl> <dbl>
## 1 1 1 combined 1.92 1.84 0.576 1 4.34
## 2 1 1 combined 1.92 1.84 0.576 2 6.57
## 3 1 1 combined 1.92 1.84 0.576 3 2.37
## 4 1 2 gesture 1.92 1.84 0.349 1 5.50
## 5 1 2 gesture 1.92 1.84 0.349 2 10.4
## 6 1 2 gesture 1.92 1.84 0.349 3 4.00
## # ℹ 3 more variables: TrialNumber <dbl>, PrevAn <dbl>,
## # Effort_Change_Ratio <dbl>Exploring structure
Let’s check in plots
So we will have to work with bimodal distribution, as we either have a
decrease in effort (Ratio < 1) or increase (Ratio > 1)
#H2: Stating causal model
We now also need a new DAG. Essentially, what we said will influence CommAtt in H1, will now also influence PrevAn because they are tightly related. For instance, more extroverted people can be expected to be better guessers, therefore the similarity of the previous answer will be higher.
## { Big5, Conc, Expr, Fam, Mod, Pcn, TrNum }The adjustment set is identical to the one of H1. Note that we are here omitting arrows going from all these variables to ComAtt. However, since PrevAn affects CommAtt, and not the other way, we anyway do not need to block this path to avoid confounds. However, if we want to asses direct effect of PrevAn on EffChange, we have to add CommAtt to the model as well
Modeling: H2
H2: A higher degree of misunderstanding will require a performer to engage in more effortful correction.
Contrast coding
- Convert collumns to factors
filtered_data$CommAtt <- as.factor(filtered_data$CommAtt)
filtered_data$Modality <- as.factor(filtered_data$Modality)
filtered_data$Participant <- as.factor(filtered_data$Participant)
filtered_data$Concept <- as.factor(filtered_data$Concept)
filtered_data$TrialNumber <- as.numeric(filtered_data$TrialNumber) # Ensure TrialNumber is numeric- Contrast-coding for categorical predictors
contrasts(filtered_data$CommAtt) <- MASS::contr.sdif(2)
contrasts(filtered_data$Modality) <- contr.sum(3)/2- Centering continuous predictors
filtered_data$TrialNumber_c <- filtered_data$TrialNumber - median(range(filtered_data$TrialNumber))
filtered_data$Familiarity <- filtered_data$Familiarity - median(range(filtered_data$Familiarity))
filtered_data$Big5 <- filtered_data$Big5 - median(range(filtered_data$Big5))- Z-scoring expressibility and similarity (PrevAn)
filtered_data <-
filtered_data |>
group_by(Modality) |>
mutate(Expressibility_z = (Expressibility - mean(Expressibility))/ sd(filtered_data$Expressibility, na.rm = T)) |>
ungroup()
filtered_data <-
filtered_data |>
mutate(PrevAn_z = (PrevAn - mean(PrevAn))/ sd(filtered_data$PrevAn, na.rm = T)) |>
ungroup()Note: maybe we should center where we z-score to get better sense of the predictors (right now, for me it's still a bit uninterpretable)Model 1 - DAG
h2.m1 <- brm(Effort_Change_Ratio ~ 1 + PrevAn_z + CommAtt + Familiarity + Big5 + Expressibility_z + TrialNumber_c + Modality + (1 | Participant) + (1 | Concept),
data = filtered_data,
iter = 4000,
cores = 4)
# Add criterions for later diagnostics
h2.m1 <- add_criterion(h2.m1, criterion = c("loo", "waic"))
# Calculate also variance explained (R^2)
h2.m1_R2 <- bayes_R2(h2.m1)
# Save both as objects
saveRDS(h2.m1, here("09_Analysis_Modeling", "models", "h2.m1.rds"))
saveRDS(h2.m1_R2, here("09_Analysis_Modeling", "models", "h2.m1_R2.rds"))
beep(5)h2.m1 <- readRDS(here("09_Analysis_Modeling", "models", "h2.m1.rds"))
h2.m1_R2 <- readRDS(here("09_Analysis_Modeling", "models", "h2.m1_R2.rds"))
# Summary
summary(h2.m1)## Family: gaussian
## Links: mu = identity; sigma = identity
## Formula: Effort_Change_Ratio ~ 1 + PrevAn_z + CommAtt + Familiarity + Big5 + Expressibility_z + TrialNumber_c + Modality + (1 | Participant) + (1 | Concept)
## Data: filtered_data (Number of observations: 2556)
## Draws: 4 chains, each with iter = 4000; warmup = 2000; thin = 1;
## total post-warmup draws = 8000
##
## Multilevel Hyperparameters:
## ~Concept (Number of levels: 21)
## Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
## sd(Intercept) 0.00 0.00 0.00 0.01 1.00 3361 3453
##
## ~Participant (Number of levels: 120)
## Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
## sd(Intercept) 0.00 0.00 0.00 0.01 1.00 3037 3838
##
## Regression Coefficients:
## Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
## Intercept 1.11 0.00 1.10 1.11 1.00 12420 6655
## PrevAn_z -0.16 0.00 -0.16 -0.15 1.00 16015 5424
## CommAtt2M1 -1.54 0.00 -1.55 -1.53 1.00 21480 4412
## Familiarity 0.00 0.00 -0.00 0.01 1.00 18103 5539
## Big5 -0.00 0.00 -0.01 0.01 1.00 17827 6080
## Expressibility_z 0.00 0.00 -0.00 0.01 1.00 15822 6089
## TrialNumber_c 0.00 0.00 -0.00 0.00 1.00 7961 6172
## Modality1 -0.00 0.00 -0.01 0.01 1.00 14413 6436
## Modality2 -0.01 0.01 -0.02 0.00 1.00 13936 6092
##
## Further Distributional Parameters:
## Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
## sigma 0.09 0.00 0.09 0.09 1.00 15108 5752
##
## Draws were sampled using sampling(NUTS). For each parameter, Bulk_ESS
## and Tail_ESS are effective sample size measures, and Rhat is the potential
## scale reduction factor on split chains (at convergence, Rhat = 1).# coefficients seem quite conservative, even for PrevAn (but remember, it's z-scored)plot(h2.m1)
# all looks ok
plot(conditional_effects(h2.m1), points = TRUE)
# the main effect seems ok, the rest of the predictors does not show nothing (maybe because of the transformation we lost the relations between them and the effort)
pp_check(h2.m1, type = "dens_overlay")
# this looks quite okay-ish
pp_check(h2.m1, type = "error_scatter_avg")
# so there seems to be quite high residual error for some extreme values (esp for the second mode)
h2.m1_R2## Estimate Est.Error Q2.5 Q97.5
## R2 0.9860988 7.752928e-05 0.9859373 0.9862448# explained variance around 98%Model 2 - log-normal
We will once again try lognormal distribution, for now leaving priors to default
h2.m2 <- brm(Effort_Change_Ratio ~ 1 + PrevAn_z + CommAtt + Modality + Big5 + Familiarity + Expressibility_z + (1 + PrevAn || Participant) + (1 + PrevAn || Concept) + (1 | TrialNumber_c),
data = filtered_data,
family = lognormal(),
iter = 4000,
cores = 4)
# Add criterions for later diagnostics
h2.m2 <- add_criterion(h2.m2, criterion = c("loo", "waic"))
# Calculate also variance explained (R^2)
h2.m2_R2 <- bayes_R2(h2.m2)
# Save both as objects
saveRDS(h2.m2, here("09_Analysis_Modeling", "models", "h2.m2.rds"))
saveRDS(h2.m2_R2, here("09_Analysis_Modeling", "models", "h2.m2_R2.rds"))
beep(5)h2.m2 <- readRDS(here("09_Analysis_Modeling", "models", "h2.m2.rds"))
h2.m2_R2 <- readRDS(here("09_Analysis_Modeling", "models", "h2.m2_R2.rds"))
# Summary
summary(h2.m2)## Family: lognormal
## Links: mu = identity; sigma = identity
## Formula: Effort_Change_Ratio ~ 1 + PrevAn_z + CommAtt + Modality + Big5 + Familiarity + Expressibility_z + (1 + PrevAn || Participant) + (1 + PrevAn || Concept) + (1 | TrialNumber_c)
## Data: filtered_data (Number of observations: 2556)
## Draws: 4 chains, each with iter = 4000; warmup = 2000; thin = 1;
## total post-warmup draws = 8000
##
## Multilevel Hyperparameters:
## ~Concept (Number of levels: 21)
## Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
## sd(Intercept) 0.00 0.00 0.00 0.01 1.00 4108 4460
## sd(PrevAn) 0.00 0.00 0.00 0.01 1.00 5318 4098
##
## ~Participant (Number of levels: 120)
## Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
## sd(Intercept) 0.00 0.00 0.00 0.01 1.00 2293 4226
## sd(PrevAn) 0.00 0.00 0.00 0.01 1.00 3211 4384
##
## ~TrialNumber_c (Number of levels: 49)
## Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
## sd(Intercept) 0.00 0.00 0.00 0.01 1.00 4800 4380
##
## Regression Coefficients:
## Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
## Intercept -0.24 0.00 -0.24 -0.23 1.00 10195 6482
## PrevAn_z -0.12 0.00 -0.12 -0.12 1.00 7996 5075
## CommAtt2M1 -1.72 0.00 -1.73 -1.72 1.00 18165 5973
## Modality1 -0.00 0.00 -0.01 0.01 1.00 15611 5899
## Modality2 -0.00 0.00 -0.01 0.01 1.00 14981 6417
## Big5 -0.00 0.00 -0.01 0.00 1.00 13219 6273
## Familiarity -0.00 0.00 -0.01 0.00 1.00 13664 6321
## Expressibility_z 0.00 0.00 -0.00 0.00 1.00 9934 6219
##
## Further Distributional Parameters:
## Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
## sigma 0.07 0.00 0.07 0.07 1.00 17435 5346
##
## Draws were sampled using sampling(NUTS). For each parameter, Bulk_ESS
## and Tail_ESS are effective sample size measures, and Rhat is the potential
## scale reduction factor on split chains (at convergence, Rhat = 1).Converting to original scale
# Extract posterior samples
posterior_samples <- as_draws_df(h2.m2)
alpha_samples <- posterior_samples$b_Intercept
# Create a list to store effects for each fixed factor
effect_list <- list()
# Helper function to calculate summary statistics
get_effect_summary <- function(effect_samples) {
mean_effect <- mean(effect_samples)
se_effect <- sd(effect_samples)
ci_effect <- quantile(effect_samples, c(0.025, 0.975))
post_prob <- mean(effect_samples > 0)
c(mean = mean_effect,
se = se_effect,
lower_ci = ci_effect[1],
upper_ci = ci_effect[2],
post_prob = post_prob)
}
# Compute expected values on the log scale
mu_1 <- alpha_samples # CommAtt 1
# Transform to original scale
effect_1 <- exp(mu_1)
# PrevAn (continuous)
if ("b_PrevAn_z" %in% colnames(posterior_samples)) {
beta_prevAn <- posterior_samples$b_PrevAn_z
effect_list$PrevAn_z <- get_effect_summary(exp(alpha_samples + beta_prevAn) - exp(alpha_samples))
}
effect_list$Intercept_mean <- mean(effect_1)
# Convert to a nicely formatted data frame
effect_df <- do.call(rbind, effect_list)
# View effects
effect_df## mean se lower_ci.2.5% upper_ci.97.5% post_prob
## PrevAn_z -0.08764297 0.0009697511 -0.08953147 -0.0857059 0.0000000
## Intercept_mean 0.78830407 0.7883040707 0.78830407 0.7883041 0.7883041Still difficult to interpret as PrevAn is z-scored
plot(h2.m2)
# looks good
plot(conditional_effects(h2.m2), points = TRUE)
# main predictor looks ok, but the rest still moderated
pp_check(h2.m2, type = "dens_overlay")
# this looks very nice
pp_check(h2.m2, type = "error_scatter_avg")
# so higher values are kind blobby with no particular trend, but the lower values are still quite errorneous
h2.m2_R2## Estimate Est.Error Q2.5 Q97.5
## R2 0.9985866 0.0001091032 0.9983392 0.9987595# explained variance 99%There are some values we should be cautious about - the relationship between predicted values and residual error is still quite correlated for lower values, and the explained variance would in real-data case probably hint to over-fitting. However, we created synthetic data that does not contain as much noise as we can expect from the real data, so it is quite likely that if we simply used all predictors we know to cause a variance here, we might be explaining most of the variation. It is quite likely this will not be the case for real data.
Splines 1 - fitting b-splines
At this point, we could be essentially content with the last model. However, as both outcome variable and the main predictor are continuous variable, we might expect nonlinear patterns in their relationship. For instance, it might be the case that for answers that are completely off, people do not feel motivated to put in more effort, and sum kind of sensible ‘threshold’ needs to be reached for the effort to increase to resolve misunderstanding, and then again decrease when the answer is very similar.
We will use a combinaiton of B-splines and Bayesian GAMs to prepare some models and try them fit to our data - even if we currently should not expect any non-linearity as we have not coded any, as we can see in the following plot
What b-splines do is they build up wiggly functions from more simple
components which are called ‘basis functions’ (B). What they do is they
divide a full range of predictor into parts and assign a parameter to
each part. The parameters are gradually turned on and off in a way that
makes their sum into a curve. Unlike polynomial regression, b-splines do
not transform predictor, but they invent a series of new synthetic
predictor variables. Each of them then exists only to gradually turn a
specific parameter on and off within a specific range of the real
predictor variable. Each of these synthetic variables is called a basis
function B. See more in Richard McElreath’s Statistical Rethinking.
We are using the code from Solomon Kurz’s adaptation of this book (https://bookdown.org/content/4857/geocentric-models.html)
# Get rid of NAs in the predictor
d <-
final_data_2 %>%
drop_na(PrevAn)
# And convert all that is necessary to factor/numerical
d$CommAtt <- as.factor(d$CommAtt)
d$Modality <- as.factor(d$Modality)
d$Participant <- as.factor(d$Participant)
d$Concept <- as.factor(d$Concept)
d$TrialNumber <- as.numeric(d$TrialNumber) Here we can see summary of each predictor
d %>%
select_if(is.numeric) %>% # Select only numeric columns
pivot_longer(cols = everything(), names_to = "key", values_to = "value") %>%
group_by(key) %>%
summarise(
mean = mean(value, na.rm = TRUE),
sd = sd(value, na.rm = TRUE),
ll = quantile(value, probs = 0.055, na.rm = TRUE),
ul = quantile(value, probs = 0.945, na.rm = TRUE)
) %>%
mutate(across(where(is.double), round, digits = 2))## # A tibble: 7 × 5
## key mean sd ll ul
## <chr> <dbl> <dbl> <dbl> <dbl>
## 1 Big5 0.98 0.57 0.09 1.93
## 2 Eff 5.1 2.79 1.43 9.79
## 3 Effort_Change_Ratio 1.35 0.75 0.28 2.19
## 4 Expressibility 0.47 0.38 0.07 1.32
## 5 Familiarity 1.13 0.55 0.15 1.9
## 6 PrevAn 0.5 0.29 0.05 0.94
## 7 TrialNumber 22.6 12.4 3 42print(d)## # A tibble: 2,556 × 11
## Participant Concept Modality Big5 Familiarity Expressibility CommAtt Eff
## <fct> <fct> <fct> <dbl> <dbl> <dbl> <fct> <dbl>
## 1 1 1 combined 1.92 1.84 0.576 2 6.57
## 2 1 1 combined 1.92 1.84 0.576 3 2.37
## 3 1 2 gesture 1.92 1.84 0.349 2 10.4
## 4 1 2 gesture 1.92 1.84 0.349 3 4.00
## 5 1 3 vocal 1.92 1.84 0.340 2 12.9
## 6 1 3 vocal 1.92 1.84 0.340 3 3.59
## 7 1 4 combined 1.92 1.84 0.719 2 6.41
## 8 1 6 gesture 1.92 1.84 0.00419 2 12.0
## 9 1 6 gesture 1.92 1.84 0.00419 3 3.44
## 10 1 8 vocal 1.92 1.84 0.591 2 9.11
## # ℹ 2,546 more rows
## # ℹ 3 more variables: TrialNumber <dbl>, PrevAn <dbl>,
## # Effort_Change_Ratio <dbl>Now we need to specify knots that function as pivots for number of different basis functions. The B variable then tells you which knot you are close to.
num_knots <- 7
knot_list <- quantile(d$PrevAn, probs = seq(from = 0, to = 1, length.out = num_knots))
knot_list## 0% 16.66667% 33.33333% 50% 66.66667% 83.33333%
## 0.0001020494 0.1540585973 0.3319149296 0.5038721841 0.6719870094 0.8383128390
## 100%
## 0.9989464423Now we can see how we chopped the data by the knots
Now we need to specify the polynomial degree which determines how
parameters interact to produce the spline.
For degree 1, two basis functions combine at each point. For degree 2, three functions combine at each point. For degree 3, four combine.
library(splines)
B <- bs(d$PrevAn,
knots = knot_list[-c(1, num_knots)],
degree = 3, # cubic spline
intercept = TRUE)This is how cubic spline with 7 knots look like
We will now add this B matrix into the data to be able to further use it
in the models
d2 <-
d %>%
mutate(B = B)
# take a look at the structure of `d3
d2 %>% glimpse()## Rows: 2,556
## Columns: 12
## $ Participant <fct> 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1,…
## $ Concept <fct> 1, 1, 2, 2, 3, 3, 4, 6, 6, 8, 8, 9, 9, 14, 15, 15,…
## $ Modality <fct> combined, combined, gesture, gesture, vocal, vocal…
## $ Big5 <dbl> 1.919899, 1.919899, 1.919899, 1.919899, 1.919899, …
## $ Familiarity <dbl> 1.8439059, 1.8439059, 1.8439059, 1.8439059, 1.8439…
## $ Expressibility <dbl> 0.576172309, 0.576172309, 0.349378655, 0.349378655…
## $ CommAtt <fct> 2, 3, 2, 3, 2, 3, 2, 2, 3, 2, 3, 2, 3, 2, 2, 3, 2,…
## $ Eff <dbl> 6.566250, 2.367010, 10.367605, 3.998448, 12.934986…
## $ TrialNumber <dbl> 2, 3, 5, 6, 8, 9, 11, 14, 15, 18, 19, 21, 22, 28, …
## $ PrevAn <dbl> 0.98450637, 0.82035722, 0.48565959, 0.09090184, 0.…
## $ Effort_Change_Ratio <dbl> 1.5116202, 0.3604812, 1.8857553, 0.3856675, 2.1021…
## $ B <bs[,9]> <bs[26 x 9]>The B matrix is now a matrix column which contains the same number of rows as the others, but has also 9 columns within that columns. Each of them correspond to one synthetic B variable.
We can now use this matrix to fit our model
h2.s1 <-
brm(data = d2,
family = gaussian,
Effort_Change_Ratio ~ 1 + B,
prior = c(prior(normal(100, 10), class = Intercept),
prior(normal(0, 10), class = b),
prior(exponential(1), class = sigma)),
iter = 4000, warmup = 2000, chains = 4, cores = 4,
seed = 4)
# Add criterions for later diagnostics
h2.s1 <- add_criterion(h2.s1, criterion = c("loo", "waic"))
# Calculate also variance explained (R^2)
h2.s1_R2 <- bayes_R2(h2.s1)
# Save both as objects
saveRDS(h2.s1, here("09_Analysis_Modeling", "models", "h2.s1.rds"))
saveRDS(h2.s1_R2, here("09_Analysis_Modeling", "models", "h2.s1_R2.rds"))Model summary
h2.s1 <- readRDS(here("09_Analysis_Modeling", "models", "h2.s1.rds"))
h2.s1_R2 <- readRDS(here("09_Analysis_Modeling", "models", "h2.s1_R2.rds"))
# Summary
print(h2.s1)## Family: gaussian
## Links: mu = identity; sigma = identity
## Formula: Effort_Change_Ratio ~ 1 + B
## Data: d2 (Number of observations: 2556)
## Draws: 4 chains, each with iter = 4000; warmup = 2000; thin = 1;
## total post-warmup draws = 8000
##
## Regression Coefficients:
## Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
## Intercept 1.35 3.28 -5.11 7.94 1.00 954 1238
## B1 0.12 3.28 -6.48 6.57 1.00 956 1246
## B2 0.16 3.29 -6.46 6.65 1.00 956 1255
## B3 0.38 3.28 -6.24 6.84 1.00 954 1248
## B4 -0.02 3.29 -6.63 6.46 1.00 956 1242
## B5 0.08 3.28 -6.53 6.54 1.00 956 1260
## B6 -0.17 3.29 -6.75 6.34 1.00 954 1242
## B7 -0.15 3.29 -6.76 6.31 1.00 958 1277
## B8 -0.21 3.29 -6.83 6.28 1.00 953 1275
## B9 -0.29 3.29 -6.92 6.19 1.00 955 1242
##
## Further Distributional Parameters:
## Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
## sigma 0.73 0.01 0.71 0.75 1.00 3635 3204
##
## Draws were sampled using sampling(NUTS). For each parameter, Bulk_ESS
## and Tail_ESS are effective sample size measures, and Rhat is the potential
## scale reduction factor on split chains (at convergence, Rhat = 1).It’s a bit difficult to see what is going on so let’s just plot it
Now with the predictor
So as expected, we see quite linear decrease similar to our previous
models
Now let’s check some diagnostics
plot(h2.s1)
# looks good
pp_check(h2.s1, type = "dens_overlay")
# because we used only main predictor, we can see quite a bad result of the posterior predictive distribution, but we will work on that
pp_check(h2.s1, type = "error_scatter_avg")
# high residual error for both low and high values
h2.s1_R2## Estimate Est.Error Q2.5 Q97.5
## R2 0.04744733 0.007830488 0.03295914 0.06351053# explained variance 4%Splines 2 - Bayesian GAMs
Now we will use a different workflow, using smooth functions with brms. brms allow for non-linear models by borrowing functions from mgcv package (REF: Wood)
This is how priors should look like for one of these functions
get_prior(data = d2,
family = gaussian,
Effort_Change_Ratio ~ 1 + s(PrevAn))## prior class coef group resp dpar nlpar lb ub
## (flat) b
## (flat) b sPrevAn_1
## student_t(3, 1.7, 2.5) Intercept
## student_t(3, 0, 2.5) sds 0
## student_t(3, 0, 2.5) sds s(PrevAn) 0
## student_t(3, 0, 2.5) sigma 0
## source
## default
## (vectorized)
## default
## default
## (vectorized)
## defaultAt this point, I will just keep priors at default
h2.s2 <-
brm(data = d2,
family = gaussian,
Effort_Change_Ratio ~ 1 + s(PrevAn, bs = "bs", k = 19),
iter = 4000,
warmup = 2000,
chains = 4,
cores = 4,
seed = 4,
control = list(adapt_delta = .99)
)
# Add criterions for later diagnostics
h2.s2 <- add_criterion(h2.s2, criterion = c("loo", "waic"))
# Calculate also variance explained (R^2)
h2.s2_R2 <- bayes_R2(h2.s2)
# Save both as objects
saveRDS(h2.s2, here("09_Analysis_Modeling", "models", "h2.s2.rds"))
saveRDS(h2.s2_R2, here("09_Analysis_Modeling", "models", "h2.s2_R2.rds"))Now we can proceed as we usually do
h2.s2 <- readRDS(here("09_Analysis_Modeling", "models", "h2.s2.rds"))
h2.s2_R2 <- readRDS(here("09_Analysis_Modeling", "models", "h2.s2_R2.rds"))
# Summary
summary(h2.s2)## Family: gaussian
## Links: mu = identity; sigma = identity
## Formula: Effort_Change_Ratio ~ 1 + s(PrevAn, bs = "bs", k = 19)
## Data: d2 (Number of observations: 2556)
## Draws: 4 chains, each with iter = 4000; warmup = 2000; thin = 1;
## total post-warmup draws = 8000
##
## Smoothing Spline Hyperparameters:
## Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
## sds(sPrevAn_1) 0.01 0.01 0.00 0.03 1.00 1526 2979
##
## Regression Coefficients:
## Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
## Intercept 1.35 0.01 1.33 1.38 1.00 10258 6031
## sPrevAn_1 -0.08 0.01 -0.09 -0.06 1.00 6284 5193
##
## Further Distributional Parameters:
## Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
## sigma 0.73 0.01 0.71 0.75 1.00 11605 5819
##
## Draws were sampled using sampling(NUTS). For each parameter, Bulk_ESS
## and Tail_ESS are effective sample size measures, and Rhat is the potential
## scale reduction factor on split chains (at convergence, Rhat = 1).plot(h2.s2)
# looks good
plot(conditional_effects(h2.s2), points = TRUE)
# again, we see similar (non-linear) trend
pp_check(h2.s2, type = "dens_overlay")
# same as before, ignoring bimodality
pp_check(h2.s2, type = "error_scatter_avg")
# still weird
h2.s2_R2## Estimate Est.Error Q2.5 Q97.5
## R2 0.04336834 0.007618562 0.02908593 0.05921277# explained variance 4%Splines 3 - GAMs with f+r effects
Now let’s use the same model, but adding our other predictors
h2.s3 <-
brm(
data = filtered_data,
family = gaussian,
Effort_Change_Ratio ~ 1 +
s(PrevAn_z, bs = "bs", k = 19) + # Smooth for Previous Answer similarity
+ CommAtt + Modality + Big5 + Familiarity + Expressibility_z + # Fixed effects
(1 | Participant) + (1 | Concept), # Random effects
iter = 4000,
warmup = 2000,
chains = 4,
cores = 4,
seed = 4,
control = list(adapt_delta = .99)
)
# Add criterions for later diagnostics
h2.s3 <- add_criterion(h2.s3, criterion = c("loo", "waic"))
# Calculate also variance explained (R^2)
h2.s3_R2 <- bayes_R2(h2.s3)
# Save both as objects
saveRDS(h2.s3, here("09_Analysis_Modeling", "models", "h2.s3.rds"))
saveRDS(h2.s3_R2, here("09_Analysis_Modeling", "models", "h2.s3_R2.rds"))h2.s3 <- readRDS(here("09_Analysis_Modeling", "models", "h2.s3.rds"))
h2.s3_R2 <- readRDS(here("09_Analysis_Modeling", "models", "h2.s3_R2.rds"))
# Summary
summary(h2.s3)## Family: gaussian
## Links: mu = identity; sigma = identity
## Formula: Effort_Change_Ratio ~ 1 + s(PrevAn_z, bs = "bs", k = 19) + +CommAtt + Modality + Big5 + Familiarity + Expressibility_z + (1 | Participant) + (1 | Concept)
## Data: filtered_data (Number of observations: 2556)
## Draws: 4 chains, each with iter = 4000; warmup = 2000; thin = 1;
## total post-warmup draws = 8000
##
## Smoothing Spline Hyperparameters:
## Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
## sds(sPrevAn_z_1) 0.00 0.00 0.00 0.00 1.00 2278 3567
##
## Multilevel Hyperparameters:
## ~Concept (Number of levels: 21)
## Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
## sd(Intercept) 0.00 0.00 0.00 0.01 1.00 3314 3910
##
## ~Participant (Number of levels: 120)
## Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
## sd(Intercept) 0.00 0.00 0.00 0.01 1.00 3210 4292
##
## Regression Coefficients:
## Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
## Intercept 1.11 0.00 1.10 1.11 1.00 11236 6961
## CommAtt2M1 -1.54 0.00 -1.55 -1.53 1.00 21880 5650
## Modality1 0.00 0.00 -0.01 0.01 1.00 15488 6225
## Modality2 -0.01 0.01 -0.02 0.00 1.00 15124 6351
## Big5 -0.00 0.00 -0.01 0.01 1.00 19293 6143
## Familiarity 0.00 0.00 -0.00 0.01 1.00 20031 5691
## Expressibility_z 0.00 0.00 -0.00 0.01 1.00 14573 5746
## sPrevAn_z_1 -0.08 0.00 -0.08 -0.08 1.00 7024 6935
##
## Further Distributional Parameters:
## Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
## sigma 0.09 0.00 0.09 0.09 1.00 21128 5955
##
## Draws were sampled using sampling(NUTS). For each parameter, Bulk_ESS
## and Tail_ESS are effective sample size measures, and Rhat is the potential
## scale reduction factor on split chains (at convergence, Rhat = 1).plot(h2.s3)
# looks ok
plot(conditional_effects(h2.s3), points = TRUE)
# now the ffect of PrevAn looks a bit different because we are back to using z-scored version
pp_check(h2.s3, type = "dens_overlay")
# looks ok but the ppd undershoots the low values and overshoot the high values
pp_check(h2.s3, type = "error_scatter_avg")
# still kind of weird, and high values positively correlated
h2.s3_R2## Estimate Est.Error Q2.5 Q97.5
## R2 0.9861186 7.755798e-05 0.9859575 0.9862589# 98% varianceSplines 4 - GAMs with lognormal distribution
As we see in the ppcheck, we are still a bit struggling with the distribution, so let’s try once again a lognormal one
h2.s4 <-
brm(
Effort_Change_Ratio ~ 1 +
s(PrevAn_z, bs = "bs", k = 19) + # Nonlinear effect for Previous Answer similarity
+ CommAtt + Modality + Big5 + Familiarity + Expressibility_z +
(1 + PrevAn_z || Participant) + # Random slopes and intercepts for Participant
(1 + PrevAn_z || Concept) + # Random slopes and intercepts for Concept
(1 | TrialNumber_c), # Random intercept for Trial Number
data = filtered_data,
family = lognormal(),
iter = 4000,
warmup = 2000,
chains = 4,
cores = 4,
seed = 4,
control = list(adapt_delta = .999,
max_treedepth = 12)
)
# Add criterions for later diagnostics
h2.s4 <- add_criterion(h2.s4, criterion = c("loo", "waic"))
# Calculate also variance explained (R^2)
h2.s4_R2 <- bayes_R2(h2.s4)
# Save both as objects
saveRDS(h2.s4, here("09_Analysis_Modeling", "models", "h2.s4.rds"))
saveRDS(h2.s4_R2, here("09_Analysis_Modeling", "models", "h2.s4_R2.rds"))h2.s4 <- readRDS(here("09_Analysis_Modeling", "models", "h2.s4.rds"))
h2.s4_R2 <- readRDS(here("09_Analysis_Modeling", "models", "h2.s4_R2.rds"))
# Summary
summary(h2.s4)## Family: lognormal
## Links: mu = identity; sigma = identity
## Formula: Effort_Change_Ratio ~ 1 + s(PrevAn_z, bs = "bs", k = 19) + +CommAtt + Modality + Big5 + Familiarity + Expressibility_z + (1 + PrevAn_z || Participant) + (1 + PrevAn_z || Concept) + (1 | TrialNumber_c)
## Data: filtered_data (Number of observations: 2556)
## Draws: 4 chains, each with iter = 4000; warmup = 2000; thin = 1;
## total post-warmup draws = 8000
##
## Smoothing Spline Hyperparameters:
## Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
## sds(sPrevAn_z_1) 0.00 0.00 0.00 0.00 1.00 2682 4039
##
## Multilevel Hyperparameters:
## ~Concept (Number of levels: 21)
## Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
## sd(Intercept) 0.00 0.00 0.00 0.01 1.00 3801 4084
## sd(PrevAn_z) 0.00 0.00 0.00 0.01 1.00 3743 3720
##
## ~Participant (Number of levels: 120)
## Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
## sd(Intercept) 0.00 0.00 0.00 0.01 1.00 1880 3012
## sd(PrevAn_z) 0.00 0.00 0.00 0.01 1.00 3870 4554
##
## ~TrialNumber_c (Number of levels: 49)
## Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
## sd(Intercept) 0.00 0.00 0.00 0.01 1.00 4249 3836
##
## Regression Coefficients:
## Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
## Intercept -0.24 0.00 -0.24 -0.23 1.00 10336 7006
## CommAtt2M1 -1.72 0.00 -1.73 -1.72 1.00 16906 6230
## Modality1 -0.00 0.00 -0.01 0.01 1.00 15101 6177
## Modality2 -0.00 0.00 -0.01 0.01 1.00 14856 6513
## Big5 -0.00 0.00 -0.01 0.00 1.00 13114 5925
## Familiarity -0.00 0.00 -0.01 0.00 1.00 13584 6054
## Expressibility_z 0.00 0.00 -0.00 0.00 1.00 9683 5907
## sPrevAn_z_1 -0.06 0.00 -0.06 -0.06 1.00 9083 7104
##
## Further Distributional Parameters:
## Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
## sigma 0.07 0.00 0.07 0.07 1.00 14570 5307
##
## Draws were sampled using sampling(NUTS). For each parameter, Bulk_ESS
## and Tail_ESS are effective sample size measures, and Rhat is the potential
## scale reduction factor on split chains (at convergence, Rhat = 1).plot(h2.s4)
# looks ok
plot(conditional_effects(h2.s4), points = TRUE)
# for main predictor ok, CommAtt has a strong affect as can be expected. The rest is probably moderated but that's okay for now
pp_check(h2.s4, type = "dens_overlay")
# but this looks very good
pp_check(h2.s4, type = "error_scatter_avg")
# so for the higher-values mode it looks ok (blob around 0), but the low values are still quite error-prone. At the same time, the scale of x-axis has quite small range
h2.s4_R2## Estimate Est.Error Q2.5 Q97.5
## R2 0.9986769 0.0001261677 0.9983885 0.9988757# 99% explained variance - we are quite likely overfitting here, but as mentioned above, this is perhaps expectable for synthetic dataDiagnostics I
Before we proceed further, let’s do first round of diagnostics
Rhat
# Extract R-hat values for each model
rhat_list <- lapply(model_list, function(model) {
rhat_values <- rhat(model)
data.frame(model = deparse(substitute(model)),
max_rhat = max(rhat_values),
min_rhat = min(rhat_values))
})
# Combine and inspect
do.call(rbind, rhat_list)## model max_rhat min_rhat
## 1 X[[i]] 1.003347 0.9999272
## 2 X[[i]] 1.003156 0.9997275
## 3 X[[i]] 1.001528 1.0003296
## 4 X[[i]] 1.002482 1.0002600
## 5 X[[i]] 1.003105 0.9996580
## 6 X[[i]] 1.004370 0.9996854All Rhat values look good
ESS
Effective sample size tells how many independent samples the model has effectively drawn from the PD. Low ESS suggests autocorrelation (i.e., sample explores one part of posterior), while high ESS means good mix
# Extract n_eff values for each model
neff_ratio_list <- lapply(model_list, function(model) {
neff_values <- neff_ratio(model) # Here we calculate ratio (not the raw number of effective samples)
data.frame(model = deparse(substitute(model)),
min_neff = min(neff_values),
max_neff = max(neff_values),
mean_neff = mean(neff_values))
})
# Combine and inspect
do.call(rbind, neff_ratio_list)## model min_neff max_neff mean_neff
## 1 X[[i]] 0.3064723 0.9110065 0.7814734
## 2 X[[i]] 0.2718387 0.9254677 0.7632834
## 3 X[[i]] 0.1190914 0.4004991 0.1717864
## 4 X[[i]] 0.1907208 0.8094037 0.5994023
## 5 X[[i]] 0.2629642 0.8886650 0.7698058
## 6 X[[i]] 0.2350103 0.9275640 0.7720078They all look good expect the very first non-linear model h2.s1 which is quite expectable since we used only main predictor
effective_sample(h2.m1) ## Parameter ESS
## 1 b_Intercept 12515
## 2 b_PrevAn_z 15913
## 3 b_CommAtt2M1 21729
## 4 b_Familiarity 17809
## 5 b_Big5 17830
## 6 b_Expressibility_z 15791
## 7 b_TrialNumber_c 8061
## 8 b_Modality1 14644
## 9 b_Modality2 13812effective_sample(h2.m2) ## Parameter ESS
## 1 b_Intercept 9992
## 2 b_PrevAn_z 7925
## 3 b_CommAtt2M1 18072
## 4 b_Modality1 15601
## 5 b_Modality2 14816
## 6 b_Big5 13249
## 7 b_Familiarity 13571
## 8 b_Expressibility_z 9882effective_sample(h2.s4) # this one probably looks the nices## Parameter ESS
## 1 b_Intercept 10266
## 2 b_CommAtt2M1 17155
## 3 b_Modality1 15003
## 4 b_Modality2 14686
## 5 b_Big5 13033
## 6 b_Familiarity 13764
## 7 b_Expressibility_z 9603
## 8 bs_sPrevAn_z_1 9100effective_sample(h2.s3) ## Parameter ESS
## 1 b_Intercept 11182
## 2 b_CommAtt2M1 21782
## 3 b_Modality1 15420
## 4 b_Modality2 15244
## 5 b_Big5 19494
## 6 b_Familiarity 19777
## 7 b_Expressibility_z 14408
## 8 bs_sPrevAn_z_1 7070LOO & WAIC
Leave-one-out (loo) validation
l <- loo_compare(h2.m1, h2.m2, h2.s1, h2.s2, h2.s3, h2.s4, criterion = "loo")
print(l, simplify = F)## elpd_diff se_diff elpd_loo se_elpd_loo p_loo se_p_loo looic se_looic
## h2.s4 0.0 0.0 3127.4 40.7 42.0 1.8 -6254.8 81.4
## h2.m2 -7.2 4.4 3120.2 40.6 34.8 1.4 -6240.4 81.1
## h2.s3 -564.5 54.3 2562.9 36.0 23.0 0.7 -5125.8 72.1
## h2.m1 -564.8 54.3 2562.7 36.2 20.1 0.6 -5125.3 72.5
## h2.s2 -5966.6 46.3 -2839.2 18.6 4.8 0.1 5678.4 37.2
## h2.s1 -5968.9 46.4 -2841.5 18.8 9.3 0.2 5682.9 37.6Here we see that both lognormal models - linear and non-linear - seem to have the best diagnostics as measured by loo
Information criterion (WAIC)
w <- loo_compare(h2.m1, h2.m2, h2.s1, h2.s2, h2.s3, h2.s4, criterion = "waic")
print(w, simplify = F)## elpd_diff se_diff elpd_waic se_elpd_waic p_waic se_p_waic waic
## h2.s4 0.0 0.0 3127.5 40.7 41.9 1.8 -6255.1
## h2.m2 -7.2 4.4 3120.3 40.6 34.7 1.4 -6240.6
## h2.s3 -564.6 54.3 2563.0 36.0 23.0 0.7 -5125.9
## h2.m1 -564.9 54.3 2562.7 36.2 20.0 0.6 -5125.4
## h2.s2 -5966.7 46.3 -2839.2 18.6 4.8 0.1 5678.4
## h2.s1 -5969.0 46.4 -2841.5 18.8 9.3 0.2 5682.9
## se_waic
## h2.s4 81.4
## h2.m2 81.1
## h2.s3 72.1
## h2.m1 72.5
## h2.s2 37.2
## h2.s1 37.6# see Solomon Kurz
cbind(waic_diff = w[,1] * -2,
se = w[,2] * 2)## waic_diff se
## h2.s4 0.00000 0.000000
## h2.m2 14.44092 8.748315
## h2.s3 1129.18036 108.542331
## h2.m1 1129.72399 108.642706
## h2.s2 11933.49758 92.630038
## h2.s1 11938.00693 92.856272Plot the comparison
Here we see identical results
model_weights(h2.m1, h2.m2, h2.s1, h2.s2, h2.s3, h2.s4, weights = "waic") %>%
round(digits = 2)## h2.m1 h2.m2 h2.s1 h2.s2 h2.s3 h2.s4
## 0 0 0 0 0 1For these data, GAMs are probably adding unnecessary complexity as we are not gaining much more explanatory power. However, we leave it open whether linear or non-linear models will be better suited for the real data, following this pipeline.
Adding priors
Similarly to H1, we will also set mildly informative priors. Let’s check what priors have been selected as default for h2.s4
# Print priors
prior_summary(h2.s4)## prior class coef group
## (flat) b
## (flat) b Big5
## (flat) b CommAtt2M1
## (flat) b Expressibility_z
## (flat) b Familiarity
## (flat) b Modality1
## (flat) b Modality2
## (flat) b sPrevAn_z_1
## student_t(3, 0.5, 2.5) Intercept
## student_t(3, 0, 2.5) sd
## student_t(3, 0, 2.5) sd Concept
## student_t(3, 0, 2.5) sd Intercept Concept
## student_t(3, 0, 2.5) sd PrevAn_z Concept
## student_t(3, 0, 2.5) sd Participant
## student_t(3, 0, 2.5) sd Intercept Participant
## student_t(3, 0, 2.5) sd PrevAn_z Participant
## student_t(3, 0, 2.5) sd TrialNumber_c
## student_t(3, 0, 2.5) sd Intercept TrialNumber_c
## student_t(3, 0, 2.5) sds
## student_t(3, 0, 2.5) sds s(PrevAn_z, bs = "bs", k = 19)
## student_t(3, 0, 2.5) sigma
## resp dpar nlpar lb ub source
## default
## (vectorized)
## (vectorized)
## (vectorized)
## (vectorized)
## (vectorized)
## (vectorized)
## (vectorized)
## default
## 0 default
## 0 (vectorized)
## 0 (vectorized)
## 0 (vectorized)
## 0 (vectorized)
## 0 (vectorized)
## 0 (vectorized)
## 0 (vectorized)
## 0 (vectorized)
## 0 default
## 0 (vectorized)
## 0 defaultFor the flat priors that have been selected as default, we can re-use our previous priors from H1,. The priors could therefore look somewhat like this
priors_h2s5 <- c(
set_prior("normal(0,0.50)", class = "b", coef = "sPrevAn_z_1"),
set_prior("normal(0,0.50)", class = "b", coef = "CommAtt2M1"),
set_prior("normal(0,0.25)", class = "b", coef = "Modality1"),
set_prior("normal(0,0.25)", class = "b", coef = "Modality2"),
set_prior("normal(0,0.25)", class = "b", coef = "Big5"),
set_prior("normal(0,0.25)", class = "b", coef = "Familiarity"),
set_prior("normal(0,0.25)", class = "b", coef = "Expressibility_z")
)Adding Interactions
Because we are currently already reaching the ceiling of R^2, we are not going to add more interactions. However, we are not excluding the option of adding few interactions when using the real data. Possible interactions include:
PrevAn x Modality - similarity affects the change in effort differently (e.g., vocal might still require more effort)
PrevAn x Expressibility - similarity in relation to effort might matter only for highly expressible concepts, and low expressible concepts are difficult to express, so also difficult to exaggerate
Familiarity x PrevAn - if the guess is really bad, only very familiar people might be motivated enough to put more effort
Big x PrevAn - same like with familiarity
Similar to the workflow adopted in H1 modelling, we would fit two new models, one with priors and one with interactions, and perform another round of diagnostics to see which model seems to have the most predictive power.